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Abstract 

X-ray grating spectra provide much spectral detail from classical nova out- 
bursts. They supplied the confirmation of continued mass loss from the nova 
in the late super-soft source (SSS) stage of the outburst. 

It is not clear a priori, what the influence of the mass loss on the spectrum 
is, apart from causing systematic blue shifts in the absorption lines. In order 
to answer this question, and to test whether it is safe to neglect this aspect 
of expansion in model atmospheres for novae in the SSS stage, physically 
consistent models for expanding nova atmospheres have been developed in 
this work. 

The very high temperatures of these models combinded with high ex- 
pansion velocities and the accompanying large radial extension make nova 
in the SSS phase very interesting objects but also physically complicated 
objects to model. In this work the general purpose radiative transport code 
PHOENIX, designed to deal with expanding atmospheres, has been chosen for 
modeling X-ray novae. 

PHOENIX has been used for this type of objects before, but careful ana- 
lysis of the old results lead to a number of new methods and improvements 
to the code, being the main achievement of this work. Firstly, essential im- 
provements to the physical treatment of NLTE have been made, including: 
new opacity expressions, a new rate matrix solver, a new global iteration 
scheme, and a new temperature correction method. Secondly, a new hybrid 
hydrostatic-dynamic nova atmosphere setup has been implemented. Thirdly, 
NLTE models are treated in pure NLTE, without LTE opacities. Also, the 
models have been made faster to compute by at least a factor 10. 

With the new framework a modest amount of models, limited by com- 
putation time, have been calculated. These models show that systematic 
results are achieved from the framework for various atmospheric conditions. 
They also show, that the influence of the expanding shell on the model spec- 
trum is important and that the model spectra are sensitive to the details of 
the atmospheric structure. The nova models are compared to the 10 well- 
exposed X-ray nova grating spectra presently available: 5x V4743 Sgr 2003, 
3x RS Oph 2006, and 2x V2491 Cyg 2008. Although the models are on a 
coarse grid and have not yet been tuned to the observations they do match 
surprisingly well. 

Also, a comparison with hydrostatic models is made. The reproduction 
of the data is clearly inferior to the expanding models. But what is more 
important is that the interpretation of the data with hydrostatic models 
leads to conclusions that are opposite to those with expanding models, e.g. 
the former requires a sub-solar O abundance and the later a super-solar. 

The models give the ability to derive accurate constraints on the physical 
conditions in the deep layers of novae that are visible only in the SSS phase. 
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Zusammenfassung 

Die grating Spektren von klassischen Novae enhalten viel spektrale Information. 
Sie haben bestaetigt, dass der Massenverlust weiterhin stattfindct im spactcn SSS 
Stadium des Nova-Ausbruches. 

Man kann nicht vorhersagen, wie der Massenverlust das Spektrum beeinflusst, 
ausser dass Absorptionslinicn blauverschoben sind. Um diese Frage zu klaeren, und 
ob man den Aspekt der Expansion wirklich vernachlaessigen kann in Model Atmos- 
phaeren fuer Novae in der SSS-Phase, wurden physikalisch konsistente Modelle fuer 
expandierende Nova-Atmosphaeren entwickclt. 

Sehr hohe Temperaturen, hohe Expansionsgeschwindigkeiten und die damit ein- 
hergehende starke radielle Ausdehnung machen Novae in der SSS-Phase zu sehr 
interessanten, aber auch physikalisch schwer modellierbaren Objekten. In dieser 
Arbeit wuerde der general-purpose Strahlungstransportcode PHOENIX verwendet, 
der fuer expandierende Atmosphaeren entwickelt wurde, um X-ray Novae zu mo- 
dcllicrcn. 

Fuer diese Objekte is PHOENIX schon eher benutzt worden, aber die gruendlichc 
Analyse der alten Modelle fuehrte zu einer Anzahl neuer Methodcn und Verbes- 
scrungen am Code, die den Hauptbcstandtcil dieser Arbeit ausmachen. Die physi- 
kalische Behandlung von NLTE wurde grundlegend verbessert mit 1) neuen Opa- 
zitaeten, einem neuen Rate Matrix Solver, einem neuen globalen Iterationsschema 
und einer neuen Temperatur-Korrketur-Methode; 2) der Implementation eines neu- 
en hydrostatisch-dynamischen Nova Atmosphacrenaufbaus; 3) der Behandlung von 
Modellen in purem NLTE, ohne LTE Opazitaeten; und 4) der Beschleunigung der 
Modcllberechnung wurde zum Faktor 10. 

Mit dem neuen Code wurde eine kleine Menge an Modellen berechnet, einge- 
schraenkt durch Rechenzeit. Diese Modelle zeigen, dass der neue Code systemati- 
sche Ergebnisse liefert fuer mannigfaltige atmosphaerische Zustacnde. Auch zeigen 
sie, dass die expandierende Huelle einen wichtigen Einfluss auf das Modellspektrum 
haben, und dass das Spektrum empfindlich ist fuer die atmosphaerische Struktur. 
Die Modelle werden verglichen mit den zehn vorhandenen gut belichteten X-ray 
Nova Beobachtungen: 5x V4743 Sgr 2003, 3x RS Oph 2006, and 2x V2491 Cyg 
2008. Obwohl das Modcllgitter grob ist und noch nicht getuned wurde auf die Be- 
obachtungen, ist die Uebereinstimmung ueberraschend gut. 

Auch wurden hydrostatische Modelle verglichen. Diese reproduzieren die Daten 
deutlich schlechter als die expandierenden. Was aber noch wichtiger ist, ist dass 
die Interpretation der Daten mit hydrostatischen Modellen zu Schlussfolgerungen 
fuchrt, die denen mit expandierenden Modellen widcrsprcchcn. Zum Bcispicl findct 
man, dass diese eine sub-solare O Haeufigkcit erfordern und jene eine super-solare. 

Die Modelle ermoeglichen es, praezise Einschraenkungen abzuleiten ueber die 
physikalischen Bedingungen in den tiefen Schichten einer Nova, die sich nur zeigen 
in der SSS-Phase. 
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Chapter 1 

Introduction 



The phenomenon of novae has long been known, the name "nova stella" 
being given by Tycho Brahe (1572). The basic principles responsible for the 
observed behavior of novae have firstly been clearly formulated in [Gro37]. 
The light originates from a continuously expanding atmosphere where the 
optical properties are determined by the density and velocity fields. Classical 
reviews arc [Pay57], [McL57] and [GS78]. Numerical models based on the 
theory of a thermonuclear runaway that cause a nova explosion, developed 
by [SS87], were found to be in excellent agreement with the optical and UV 
observations available at that time. A recent collection of work on the topic 
of classical novae is [BE08]. The following short summary is based on these 
publications and the references therein. 



1.1 Novae: the classical theory 

The commonly accepted model for a nova is a close binary system with one 
member being a white dwarf and the other a larger, cooler star that fills 
its Roche lobe. Because it fills its lobe, any tendency for it to grow in size, 
because of evolutionary processes or for the lobe to shrink because of angular 
momentum losses, will cause a flow of gas through the inner Lagrangian 
point into the lobe of the white dwarf. The size of the white dwarf is small 
compared to the size of its lobe and the high angular momentum of the 
transferred material causes it to spiral into an accretion disk surrounding 
the white dwarf. Some viscous process, as yet unknown, acts to transfer mass 
inward and angular momentum outward through the disk so that a fraction 
of the material lost by the secondary ultimately ends up on the white dwarf. 
Over a long period of time, the accreted layer grows in thickness until the 
bottom reaches a temperature that is high enough to initiate thermonuclear 
fusion of hydrogen. Given the proper conditions a thermonuclear runaway 
(TNR) will occur at the bottom of the accreted hydrogen-rich envelope. 
Once the energy released by the thermonuclear reactions has reached the 
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surface of the white dwarf its effective temperature rapidly increases. With a 
luminosity of the order of the Eddington luminosity and the radius of a white 
dwarf, temperatures at peak will exceed 10 6 K, depending on the mass of the 
white dwarf [STWS96]. The nova becomes a very strong X-ray emitter with 
a soft spectral energy distribution of a hot stellar atmosphere. After ejection 
of the nova shell, which is expanding and cooling adiabatically, the temper- 
ature of the expanding effective photosphere rapidly drops. When hydrogen 
starts to recombine the X-ray flux decreases rapidly as the expanding enve- 
lope becomes opaque to X-rays within a few hours. At this stage, called the 
'fireball phase', the maximum in visual luminosity is reached. 

After the rise to visual maximum the shell continues to expand so that 
the density and opacity decrease and the effective photosphere moves in- 
ward. In the initial outburst only a part of the accreted material is ejected. 
The part of material that doesn't reach velocities higher than the escape 
velocity slows down, and gradually turns back into hydrostatic equilibrium 
on top of the white dwarf surface, providing enough fuel for hydrogen burn- 
ing during the so-called 'phase of constant bolometric luminosity'. With the 
constant luminosity and the effective photosphere contracting the effective 
temperature rises to a level where a second phase of X-ray emission begins, 
the so called 'supersoft source phase'. 

At some point in time the hydrogen burning turns off. The atmosphere 
of the white dwarf is then expected to cool at constant radius, and as a 
consequence the X-ray luminosity would drop. However, this 'turn-off phase' 
is, as yet, the least studied and least understood. A possible explanation 
could be that due to mass losses in a stellar wind the burning shell becomes 
too thin to sustain the hydrogen burning. 

Two examples of lightcurves measured in different wavelength bands for 
recent novae are shown in figures 1.1 and 1.2. 

1.2 The advent of 'high-res' X-ray spectra 

Before the advent of the 'high-resolution' X-ray spectra from the gratings on 
board XMM-Newton and Chandra theory and observations were in excel- 
lent agreement. The lightcurves could be reproduced with evolution mod- 
els based on a TNR process [SS87] and the optical and UV spectra were 
reproduced with radiative transport models [HSS + 95]. A classification of 
objects could be made based on the time scales of typical developments in 
the lightcurves of novae (the speed classes). Differences in chemical compo- 
sitions for the speed classes could be determined from the spectra and these 
classes could be related to distinct initial conditions for the TNR. 

With the very unexpected results from the grating spectra this picture 
changed. Many of the details revealed in these X-ray spectra cannot yet 
be explained by the established theory. So far, each new X-ray observation 
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20 40 60 80 100 

days after outburst (2006, Feb 12.84) 

Figure 1.1: Light curves of the classical nova RS Oph 2006 for X-ray (Swift XRT) 
and optical (AAVSO amateurs) wavelength bands. Unfortunately, the simultaneous UV 
observations made with Swift are overexposed. About 32 days after visual maximum 
the X-ray count rate increases strongly, and is initially very variable. Around day 45 it 
stabilizes soon after which it starts to decline gradually, (from [Nes09]) 

of a nova revealed new facets and no nova has shown the same behavior 
as preceding ones. Therefore, no classifications could yet be made, which 
gravely complicates the development of theoretical descriptions. 



1.3 This work 

One important facet in the development of the theory is detailed spectral 
analysis. The comparison of synthetic spectra with observational data can 
yield a lot of insight in the structure of the atmosphere that produces the 
radiation. PHOENIX is a general-purpose, state-of-the-art stellar and plane- 
tary atmosphere code [Hau92], [HB99], [HBBA03], [HB04]. It can calculate 
atmospheres and spectra of stars all across the HR-diagram including main 
sequence stars, giants, white dwarfs, stars with winds, TTauri stars, Novae, 
Supernovae, brown dwarfs and extrasolar giant planets. The universal ap- 
plicability makes the code well tested for a large diversity of atmospherical 
conditions. PHOENIX is one of the most advanced atmosphere codes that is 
capable of calculating expanding, optically thick atmospheres. 
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Figure 1.2: The simultaneous lightcurves of V2491 Cyg 2008 have a spectacular time 
resolution and accuracy for three wavelength bands. The supersoft X-ray phase lasts from 
day ~36 to day ~53. Note the steep decline in the X-ray lightcurve after its maximum 
and the difference to RS Oph 2006 in figure 1.1. (from [Nes09]) 



PHOENIX has been applied before to X-ray novae for the object V4743 
Sgr [PHNS05] . In the process of finding suitable models for other nova ob- 
servations (RS Oph 2005 and V2491 Cyg 2008) a number of new methods 
and improvements to the code have been developed. With these new de- 
velopments a framework is obtained from which it is possible to compute 
physically realistic, massive NLTE models for geometrically extended and 
expanding nova atmospheres. This has been the major achievement of this 
work. In the little time left after completion and verification of the frame- 
work, a small number of models have been computed. Although just a few 
first results were obtained from the new models, the agreement with the 
high-res X-ray observations is good. Also, large improvements in the spec- 
tral fits to observations of V4743 Sgr are obtained with respect to results 
from the old code. 

After an overview of the theory of radiation transport and "stellar" at- 
mosphere modeling (chapter 2), the new construction procedure of the nova 
atmospheres used in this work is described (chapter 3). Then a detailed 
description of the new methods developed in this work is given (chapter 4) 
and the first results obtained with the new models are shown (chapter 5). 



Chapter 2 

Atmosphere modeling 



2.1 Radiation transport basics 

The definitions given here comprise only some of the most important radia- 
tion quantities, namely those which are needed in the following sections. It 
is mainly a summary of [RL79]. Another useful description of the processes 
that form and affect the light a stellar atmosphere emits is given in [Rut95] . 

In the following, all quantities subscripted with A denote wavelength 
dependence. 

Intensity Light can be described by rays that follow geodesic paths. The 
specific intensity I\ of a ray, at wavelength A, is defined as 



with dE\ the amount of energy transported through the area dA at the 
location r, with n the normal to dA, between times t and t + dt, in the 
frequency band between A and A + dA, over the solid angle dfi = sin 6 d0d(fi 
around the direction I with spherical coordinates 6 and (p. 

Emissivity The contribution of the local emissivity rj\ (defined per cm 3 ) 
to the intensity of a ray is 



where s is the path length along the ray. There are different sources of 
emissivity, as explained in section 2.5. 

Extinction coefficient The extinction coefficient xx specifies the energy 
fraction taken from a ray per unit path length. It can be interpreted as a 
geometrical cross-section per unit volume. 



dE x = h (f, I, t) (I, n) cos 9 dA dt dA dQ, 



(2.1) 



dl\(s) = Vx(s)ds 



(2.2) 



(2.3) 



5 



6 



CHAPTER 2. ATMOSPHERE MODELING 



There are different processes that cause radiative extinction, as explained in 
section 2.5. 

The radiative transfer equation The rate of change of the intensity 
of a ray is given by the combined effects of emission (equation (2.2)) and 
extinction (equation (2.3)) 

^ = m-xxh (2.4) 

This basic equation expresses that the intensity along a ray changes if pho- 
tons are added or taken from it. Generally this monochromatic differential 
equation is complicated by the fact that the opacities (xx and r]\) depend 
in a direct and/or indirect way on / for all wavelengths. This is described 
in section 2.5 and 2.9. 

2.2 Common radiation definitions 

In the process of solving the radiative transfer equation (2.4) the following 
common definitions are important. 

Source function The source function is the ratio of the emissivity (equa- 
tion (2.2)) and extinction coefficient (equation (2.3)) 

Sx = rix/xx (2.5) 
Mean intensity The mean intensity J\ is the angle integral of I\ 

Ja = h I h dn = \ S\ h dfi (2,6) 

In the second step ID spherical coordinates are used (see equation (2.1)) 
and \i = cos 6. J is called the zeroth angular moment of the intensity. 

Flux The flux H\ is the energy flow outwards, in radial direction f 

Hx ^hJ h{r ' f)dn = \ J\ ^ h d ^ (2 - 7) 

where I is the unit vector pointing from the coordinate center into direction 
(1, (f>, 9). H is called the first angular moment of the intensity. 

K integral In analogy to the other moments, the second angular moment 
of intensity K\ is defined as 



1 



2 



Kx = ^ / ^hdfi (2.8) 



i 



K is proportional to the radiation pressure. 
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Optical length and optical depth The optical length is defined as 

dTx = X\ds (2.9) 

This definition is valid for any line of sight. A special case is the radial 
optical depth 

oo 

rx(r) = J X x(r)dr (2.10) 

r 

which measures the optical depth along the radial line of sight from t\ = 
at the observer's eye located at r = oo. 

2.3 Thermodynamic equilibrium 

2.3.1 Radiation 

In thermodynamic equilibrium (TE) the radiation field in an adiabatic en- 
closure is isotropic and its wavelength dependence is given by the Planck 
function B [Uns38], that is parameterized by the temperature T 

IF = B,(T) S ^^j^ (2.11) 

2.3.2 Matter 

Just like the radiation, the energy distribution of matter depends on the 
temperature T only. Each atomic species (chemical element) has its own 
energy-level distribution described in the following. 

Excitation 

The electronic excitation of atoms is described by the Boltzmann equation 



i^is \ 9js 

— I = exp 

n is J TE 9is 



Xjs - Xi 
kT 



9js 

exp 

9is 



hc/X 



kT 



(2.12) 



where rii s is the number density (cm -3 ) of the atom under consideration in 
ionization stage s and in excitation level i, g is the statistical weight of the 
level, x the excitation energy measured from the ground state, and hc/Xij 
the energy of the photon that corresponds to the energy difference between 
the levels. 

The total TE number density iV™ of the atom in ionization stage s is 
then the sum over all excitation levels 

TE TE 

Nj E = J>™ = ^E^ e " WfcT = ^Q^(T) (2.13) 

; 90S ; 90S 



8 



CHAPTER 2. ATMOSPHERE MODELING 



where Os denotes the ground state of ion s. The TE partition function 
Qj E (T) is defined as 



(2.14) 



Ionization 

The ionization balance is described by the Saha equation 
'N k+1 \ 2 Q k+1 f 2irm c kT\ 3/2 ^_ Xk/kT 



N k J TE N c Q k \ h 2 J 



(2.15) 



where N c is the electron number density and m e the electron mass. The 
factor of 2 is the statistical weight of the electron. If chemical element E 
has n ionization stages, then equation (2.15) gives n — 1 equations. The 
non-linear system of equations is closed by the constraint 



n e = ^2n s 



(2.16) 



Solving this system is straight forward. One starts with a guessed value 
for Aq, uses the Saha equation for the next higher N s recursively. The 
summation on the right-hand side of equation (2.16) yields the guess-related 
value, say a. Finally, all N s are scaled by the factor Ne/cl. 

The ionization balance is updated for all chemical species, which yields 
a new value for the electron density. And since the ionization balance de- 
pends on the electron density, the system of equations (2.15) and (2.16) is 
solved iteratively for each species, until the electron density is converged to 
a prescribed precision. 



Saha-Boltzmann distribution 

The combination of equations (2.12) and (2.15) gives the TE population 
density of a particular level % of ionization stage s 



n 



TE _ TE Sis N e 



= n 



9c 



h 2 



27rm c kT 



\ 3/2 

J exp 



Xis 



Xc 



kT 



(2.17) 



where cs denotes the "continuum" state to which the ionization takes place. 
In this work it is assumed that all ionization and recombination processes 
involve the ground state g of the next higher ionization stage n cs = n 9jS+ i. 
Therefore, all rii s (of the element under consideration) can be computed 
with equation (2.17), starting from the uppermost continuum state, when 
the ionization balance (equation (2.16)) is solved. The uppermost continuum 
state is the bare nucleus. 
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2.3.3 Interaction of radiation with matter 
Detailed radiative balance 

Since in TE the radiation field inside a closed system is given by the Planck 
law, equation (2.11), it must be reproduced by interactions of the radiation 
field with matter. This condition is called detailed radiative balance. It 
states that in TE the number of photons which are locally absorbed in the 
momentum range \p,p + dp\ must be equal to the number of photons emitted 
in that range. 



Detailed rate balance 

Another balance can be derived for TE, based on the Saha-Boltzmann dis- 
tribution, the detailed rate balance. It states that in TE the number of atoms 
leaving the excitation state % must be equal to the number that enter state 
i. 



KirchhofT-Planck relation 

From detailed radiative balance it follows that 

VI E = xTh (2.18) 

A classical interpretation of the interaction between radiation and matter is 
that the interaction processes can be divided into pure scattering and true 
absorption + emission processes. Here 'true' means that in the process of 
absorption the photon is destroyed and its energy is added to the thermal 
energy of the gas. In contrast, in pure scattering processes the photon inter- 
acts with the scatterer and propagates in a new direction. This interaction 
can be elastic or inelastic. In the latter case a change in the scatterer's 
internal excitation state is produced, but the assumption for scattering is 
that no energy of the photon is converted to kinetic energy of the particles 
in the gas. 

In this classical picture the macroscopic emissivity rj and extinction co- 
efficient x can be- written as 

r]x = r ] f + a x I x (2.19) 
XX = k\ + <?\ (2.20) 

where 'th' means thermal, k is the true absorption and a the pure scattering 
coefficient. 

Inserting these rj and x i n equation (2.18) the pure scattering terms 
cancel and one obtains 

^ h = K\B\ (2.21) 



This is called the Kirchhoff- Planck relation. 
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2.4 Local thermodynamic equilibrium 

In stellar atmospheres there is a net flux outwards H(r) / so that the 
assumptions of TE are not met. Local thermodynamic equilibrium (LTE) is 
a less stringent requirement than TE, because the radiation field is allowed 
to depart from the isotropic TE field (equation (2.11)), as long as it does not 
change the state of matter. This may occur for example if locally collisional 
processes play a dominant role over photon processes. 

Unlike in the TE case, the radiation field is not known a priori. Generally 

h ± B X (T) (2.22) 

but I\ is described by the radiation transport equation (2.4). 

For this reason detailed radiative balance is not valid in LTE. The 
Kirchhoff-Planck relation (2.21) does not hold exactly. However, it is a 
reasonable approximation (with B\ replaced by I\) if the gradient of the 
radiation field is small over the mean free path of a photon, which is true in 
LTE. 

In LTE the local state of matter is the same as for TE, given by equation 
(2.17), but now with a local temperature T 



„LTE _ TE 
'Hs ~ n is 



LTE 9is N e ( h 2 \ 3/2 Yio-Y^l ( 2 - 23 ) 



exp 



Xis - Xc 

kT 



cs g cs 2 \2irm c kT 
In LTE detailed rate balance still holds. 

2.5 Radiative rates and opacities 

There are different sources of opacity, with physically different underlying 
processes. The radiative processes accounted for in this work are: 

• Atomic Bound-Bound (bb) processes (excitation, deexcitation) 

• Atomic Bound-Free (bf) processes (ionization, recombination) 

• Atomic Free- Free (ff) processes (Bremsstrahlung) 

• Thomson scattering 

The emissivity (equation (2.2)) and the extinction coefficient (equation (2.3)) 
can therefore be written as 

Vx = T& h + Ttf + 4 + <TxJx (2.24) 
Xx = xf + Xx" + xl + <?x (2.25) 
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where a is the scattering coefficient. In this work only Thomson scattering 
is considered. This is an elastic scattering process. Compton and inverse- 
Compton scattering are inelastic processes that occur for high-energy pho- 
tons and high-energy electrons respectively. A self consistent treatment 
would introduce an inter-wavelength coupling of the radiation field. This 
gravely complicates the requirements for the radiative transport method and 
are therefore not treated in this work. Rayleigh scattering can be assumed 
to be negligible for hot stellar atmospheres, see [Boe89]. 

Equations (2.24) and (2.25) are the macroscopic quantities describing 
the interaction of the radiation field with matter. 

2.5.1 Thomson scattering 

In general Thomson scattering depends on the angle between the incident 
and the outgoing photon, but usually, as in this work, it is assumed to be 
isotropic. This simplification allows to write the contribution of scattering to 
the local emissivity rj\ in equation (2.24) proportional to the mean intensity 



The amount of radiation scattered out of the ray is a\I\, and a fraction 
of the local mean intensity <J\J\ is scattered into the ray. For small values 1 
of 7 = (hc/\)/(moc 2 ) <C 1 the Thomson scattering coefficient is [Pom73] 



with n e being the electron density, e the electron charge and mo the electron 
mass. 

2.5.2 Bound-Bound 

There are three forms of radiative bound-bound (bb) processes: absorption, 
induced emission and spontaneous emission. The radiative rates R at which 
these processes occur (per second and) per particle are expressed in the 
Einstein coefficients Bji and Aji, where i is the lower and j the upper 
level of the transition [MW84]. 

Radiative rates 

The absorption rate per particle is proportional to the local radiation field 



with the absorption profile <p. The profile function is area normalized 



Jx- 




(2.26) 




(2.27) 




(2.28) 



x If A = lA then 7 has the numerical value 7 = 0.024. 
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Note that B^ is defined as the rate per units of intensity in the frequency 
representation 2 . The reason for adopting this definition, and the comparison 
with a definition based on the wavelength representation of the intensity, is 
discussed in section A. 3. 

Conversion of equation (2.27) from the frequency representation to the 
wavelength representation is achieved by using the equivalence of the inte- 
grated radiation fields in both representations 

/•oo roc rO j poo 

/ J A dA= / J u dv = J,TvdA= / -^JudX (2.29) 

JO Jo Joe dA Jo X 

Since the right and the left hand side both are integrals over the same 
variable A their integrands must be equivalent. The same holds for the 
integrated profile functions in both representations 

(fi\d\ = / 4> v dv= / 4> v -,d\ (2.30) 
'o Jo Jo X 

The representation conversion rules are thus found to be 

Jx = ^Jv (2.31) 
<Px = (2.32) 

where v = c/A is substituted in the expressions for the frequency represen- 
tation. Note that the profile function is still normalized after the change of 
the integration variable. 

Using equations (2.31) and (2.32) the absorption rate equation (2.27) 
can be written in wavelength representation as 

ro ° A 2 , A 2 , du , . „ f°° , A 2 
'o 



/•OO VZ VZ Jj, /"OO yZ 

Rij = Bij / — ( j )x — J x — dX = Bij / (f>x—Jx dX (2.33) 
Jo c c a\ Jo c 



The emission rate per particle consists of a spontaneous emission term 
and an induced emission term. The latter is proportional to the local radi- 
ation field. 

/•oo \ 2 

Rji = Rfr + R l f = Aji + B jt / tpx Jx dX (2.34) 

Jo c 

Here again Bji is defined as the rate per units of intensity in the frequency 
representation, and the factor A 2 /c comes from the conversion to the wave- 
length representation. 



2 In the frequency representation the analogon of equation (2.1) is AE V = 
I v (r,l,t)(l, n) cos 6 dA dt dv dQ, 
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These rates R are total rates per particle, where total means integrated 
over all angles and all wavelengths. 



R ij = J ,% A (A,0)dOdA (2.35) 

R jt = y% iA (A,ft)dfMA (2.36) 

The angle and wavelength dependent radiative rates and &ji,\ are 
given by 



A 2 

%,A = -^B^x—h (2.37) 



%,A = + = ^ (^a + Bji^I^ (2.38) 

It can be shown from quantum mechanical considerations, that the profiles 
for induced emission and spontaneous emission must be the same [Dir47]. 

fx = V>A (2.39) 
Einstein relations 

In TE the radiation field must be Planckian. This forces the rates to satisfy 
detailed radiative balance (see section 2.3.3). The rate at which photons are 
removed from the field must be equal to the rate at which photons are added. 
In general, at frequency v all bound-bound processes, i.e. the transitions 
between all levels i and j, interact with the radiation field according to their 
profile functions 4>%j,v an d ipij,v (i n t ne following simply written as <p u and 
^v). Bound-free and other processes can be neglected for the moment, as 
justified below. The detailed radiative balance condition thus becomes 

^ riiBij4> v J u = ^ rijAjiipv + ^2 njBjiip v J u (2.40) 

ij ij ij 

Solving for J and inserting the TE radiation field J v = B v yields 

If one assumes that the profile functions are narrow, then in the line centers 
i>ij the summations over all transitions reduce to only one line. It can be 
shown [Mih78], that in TE for each transition the profile functions can be 
assumed to be ipij,v = <t>ij,w Thus it follows 

g _ Aji 1 _ Ajj 1 (2 42) 

Vii Bji n * B v _ 1 Bji SiBjj c hu ;j /kT _i 1 ' ' 

rijBji 9j Bji 
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where in the last step the Boltzmann level distribution, equation (2.12), 
is used. The equivalence to the Planck function (2.11) yields the Einstein 
relations 



These relations reflect, that the Einstein coefficients are properties of the 
atoms only and are independent of the radiation field. 

Formally, other interactions of the radiation field, like bound-free and 
free-free transitions, must be considered in the detailed radiative balance 
equation (2.40). But in general, since the derivation holds for any gas at any 
temperature, one can find a situation, where the other processes are negligi- 
ble at the exact center of a sharp bound-bound transition line. Furthermore, 
if the Einstein relations (between atomic properties, and their independence 
of the radiation field) are valid for strong bound-bound transitions, then it 
must be valid for weaker transitions too, since their interaction with the 
radiation field is not different. 

Using the Einstein relations (equations (2.43) and (2.44)) the downward 
rates can be written as functions of B, L j , the same Einstein Coefficient as the 
upwards rates 



Opacities 

The bound-bound opacities of equations (2.24) and (2.25) are in fact macro- 
scopic rates at which energy is removed or added to the ray. They can be 
expressed as the microscopic photon rates per particle multiplied by the 
particle density n and the energy per photon. However, the macroscopic 
picture in equation (2.4), is that the extinction depends on the intensity of 
the ray / and the emission does not depend on /. In order to match the 
two descriptions, the induced emission rate is treated as negative absorption 
instead of positive emission. 




(2.43) 



(2.44) 




(2.45) 



(2.46) 




(2.47) 



(2.48) 
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Inserting the expressions for the rates, equations (2.37) and (2.46), the emis- 
sivity and the extinction coefficient become 

Qi 2/ic , . . 

^ (2.49) 



„bb _ 
Vij,X - 


he 




ny 






he 




A 2 




c 



9j Aj] 

ni<f>\- rij—^x) (2.50) 
9 j J 

and the bound-bound opacities of equations (2.24) and (2.25) are given by 
the sum over all possible transitions ij. 

Using these expressions for the line opacities, the single line source func- 
tions 5 b A can be written as 

ebb _ _ 2/ic 2 / ntg^x 



(2.51) 

2hc 2 fbi -Ji£- 

1 e A fcT _ I 



~ A 2 \bj 

where in the last step complete redistribution = <p\ was assumed. When 
LTE holds, then in the line center the single line source function adopts 
the Planck value Sf^^ = B\ . In the wings, the contribution to the (total) 
source function S\ is very small, so that generally, in the continuum the 
source function is dominated by the bound- free terms. 



Line profiles 

The energy of a bound-bound transition is not sharp. It is expressed in the 
rate for such a transition with a profile function with a specific width. The 
profile function is normalized, meaning that the rate originates from a range 
of wavelengths and that the integral over the range yields the total rate. 
The precise form of the profile functions has a very important influence on 
the run of the total opacity over wavelength (equations (2.24) and (2.25)), 
and thus on the radiation field. 



Redistribution In the classical interpretation of the interaction between 
radiation and matter, described in section 2.3.3, the bound-bound transi- 
tions are divided into pure scattering and true absorption + emission pro- 
cesses. In the true absorption + emission processes there is no correlation 
between the incoming and outgoing photon. 

In the pure scattering processes both the direction and the frequency 
of a photon may change. These changes are described by a redistribution 
function R [Mih70] 

/ i?(A', A, n', n) dA' dA^-^ = 1 (2.52) 
Jo 47r 47r 
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where R(\', A, ft', ft) dX' dA^-^p gives the probability that a photon will be 
scattered - from direction n' in solid angle dW and wavelength range (A', A'+ 
dA') - into solid n in solid angle dTL and wavelength range (A, A + dA). The 
redistribution function R contains within it both a normalized absorption 
profile 4>\' and a normalized emission profile V'A The absorption profile 4>y is 
obtained by integration of R over the outgoing wavelengths and solid angle, 
and the emission profile ip\ by integration of R over the ingoing wavelengths 
and solid angle 

</> x ,(tf) = j> J i?(A',A,n',n)dA^ (2.53) 

^ A (n) = j J R(X',X,n',n)dX'^- (2.54) 

Equation (2.52) gives the full angle-frequency dependence of the emission 
profile. It is difficult to treat the radiative transfer problem in the degree of 
generality implied here, but some useful simplifications of the problem can 
be made. 

For example, one can assume a angle independent redistribution func- 
tion R(X',X). From equations (2.53) and (2.54) and the normalization by 
equation (2.52) it follows that 

R(X',X) = M X (2.55) 

In the special case that R(X' , A) = R(X, X') this reduces to 

<f> x = ipx (2.56) 

which is called complete redistribution. A good approximation to this case 
occurs, for example, when the time of interaction is long enough for collisions 
to occur. These redistribute the exited electron to the degenerate sub states 
of the upper level and thus completely remove the angle and frequency 
correlation to the initial absorption process. 

A sophistication of complete redistribution is partial redistribution, e.g. 
[Mih78]. Such methods have the disadvantage that they introduce another 
direct inter-wavelength coupling of the radiation field, which tremendously 
complicates the radiation transport problem (section 2.9). 

In the microscopic description of this section 2.5.2 the clear distinction 
between true absorption + emission and pure scattering processes cannot 
be made. Every interaction, either absorption or emission, is treated indi- 
vidually based on its profile function. So the profile functions statistically 
account for all redistribution effects. Generally, such profile functions should 
be derived from quantum mechanics. 

In this work the simplifying assumption of complete redistribution is 
made and the profile functions <p and tjj are assumed to be given by a Voigt- 
profile. 
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Voigt profile The profile functions used in this work are defined in the 
frequency representation, which is consistent with the definitions of the Ein- 
stein coefficients (see equation (2.27)). A comparison with other definitions 
is given in the appendix (section A). The Voigt profile is given by 

H(a, u) ,„ „„. 

^(a,u)=^(a,u)= ( > (2.57) 



- ^ / out/ yikT/m + 1 

vd = —^/2kT/m + i 2 = 

c Aq 



(2.58) 



a 



with the dimensionless frequency offset u and the damping parameter a 
determined by 

u = = C/X ~ C/A ° (2.60) 
Av D Au D 

a = n = — a (2-61) 

Avd being the Doppler width, £ the microturbulence parameter, and the 
dispersion width Av<i determined by the parameter T is given by 

Au d = r/(47r) (2.62) 

H is called the Voigt function. The Voigt-profile results from the convolution 
of a Gauss profile with width Avd 

i,(Au D ,u) = ^J-—e~ u2 (2.63) 



and a Lorentz profile with width Av^ 

}__Aim 

it Av 2 + Av 



The Gaussian part in the Voigt profile describes line broadening by thermal 
motions of the atoms, which causes Doppler shifts in the wavelength of ab- 
sorbed and emitted photons. The Lorentzian part is the exact profile shape 
resulting from radiative damping, which is caused by the limited lifetime of 
excited quantum states. Currently the other broadening processes, the Van 
der Waals broadening and the quadratic Stark effect, are all described by 
Lorentz profiles. This simplification allows to use the fact that the convolu- 
tion of two Lorentz profiles again yields a Lorentz profile with a width equal 
to the sum of the original widths. 

Collision broadening is a very wide field of work, quoting Rob Rutten: 
"A physicist may easily spend a whole career on collisional line broadening" 
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[Rut95]. Extensive work on line broadening by electron collision processes 
was done by [Gri74]. 

In the inner regions of the atmospheres modeled in this work, where the 
electron densities are high and the Stark effect is a very effective broadening 
process, the current description yields T-factors that are too small. In the 
outer regions the radiative damping factor gets more important, so that the 
approximations for the collisional damping are less important there. And 
for the shapes of individual lines the outer region of the atmosphere is more 
important than the inner region. However, the opacity in the inner regions 
do influence the overall structure of the atmosphere, so that the indirect 
influence of the too small T-factors might yet be important. 

2.5.3 Bound-Free 
Radiative rates 



Not yet... 
Opacities 

In analogy to the bound-bound opacities the bound-free opacities follow 
from the rates. 



Not yet... 

2.5.4 Free-Free 

Free-free transitions are transitions between unbound states. The rates of 
these transitions do not influence the population numbers of the bound 
states. Thus they do not influence the statistical equilibrium (see section 
2.8.1) and the rates do not need to be computed. 

The free-free opacities are calculated from Gaunt factors from [Sut98] 
and [ISK + 00] with routines from the publicly available CHIANTI 4 database 
(see section 2.7). These opacities are assumed to occur at TE rates, so that 
the emissivities follow from the Kirchhoff-Planck relation (2.21). 




(2.65) 
(2.66) 



Xci — ••• 



Vic = ■■■ 



(2.67) 
(2.68) 
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2.6 Collisional rates 

The collisional rates Cjj are independent of the radiation field. They are de- 
termined from theoretical collision strengths that are input data to PHOENIX. 
The data are part of the atomic databases described in section 2.7. These 
contain data for electronic collisions and protonic collisions. 

The collisional rates occur at thermodynamic equilibrium values as long 
as the gas of colliders has a Maxwellian velocity distribution [Mih78]. This 
is a reasonable assumption for the atmospheric conditions treated in this 
work. Therefore, the contribution of the collisional rates to the total rates 
Pij, equation (2.70), always drives the population numbers towards the TE 
(Saha-Boltzmann) distribution. 

2.7 Atomic data 

The calculation of the X-ray emission and the spectral features of a hot 
stellar plasma requires atomic transition rates and energy levels of highly 
ionized chemical elements. In PHOENIX multiple high energy databases are 
available 

• CHIANTI 3, see [DLM+97] 

• CHIANTI 4, see [DLYD01] 

• CHIANTI 5, see [LDY+06] 

• APED 1.3.1, see [SBLR01] 

In this work primarily the CHIANTI 5 database has been used, since it is 
the most modern and most complete of the four. For many ions the APED 
database contains only CHIANTI data. A systematic comparison of the 
different databases with PHOENIX was made by [Pet05]. 

2.8 Non-LTE 

The assumption of LTE is more realistic than TE, but is not generally valid 
in stellar atmospheres. When the assumption of LTE is dropped, it is called 
Non-LTE (NLTE). In NLTE the state of matter is no longer given by the 
Saha-Boltzmann distribution and is no longer fixed by the local temperature. 

2.8.1 Statistical equilibrium 

In this work statistical equilibrium is assumed. This means that the popu- 
lation numbers and the radiation field adjust so quickly, that the upwards 
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and downwards rates are in balance. This allows to write down the rate 
equations for a snapshot of the atmosphere 

= ^ = E n i p * - n * E p ^ ( 2 - 69 ) 

which describe the net change of occupation number n,i of level i depending 
on the transitions (bound-bound and bound-free) from and to all possible 
levels j. Pij denotes the rate of the transition from i to j per particle per 
second. The rates P contain the radiative R and the collisional rates C 

Pij = Rij + Cij (2.70) 

The radiative rates are given by the sum over all bound-bound (equations 
(2.33) and (2.45)) and bound- free rates (equations (2.65) and (2.66)). 



2.8.2 Rate matrix equation 

One rate equation (2.69) exists for each state i. Its variables are all rij with 
j / i. This is a linear system with a dimension equal to the number of levels 
N in the atom. However, the system is not yet of maximum rank because 
one equation factors out to an identity. The system is closed by replacing a 
redundant equation with the constraint of number conservation 

N 

N E = Y,n l (2.71) 

i=i 

Here the summation is performed over all the N levels of the atom under 
consideration in all ionization stages. All levels of all stages are numbered 
continuously, starting from the ground state of the neutral atom. In this 
notation the last level N is the bare nucleus. The matrix equation can be 
written in the general form 

Pn = b (2.72) 

with P being the rate matrix, n the vector of occupation numbers (number 
densities) to solve for, and b the vector that has only one non-zero (Ne) 
element from the constraint equation. An explicit rate matrix is shown 
in section 4.6, figure 4.11, for the example of a simplified atom with four 
ionization stages, the last stage being the bare nucleus. 

This simple treatment changes the ionization balances of all elements 
independently. This results in a change of the electron density. Because all 
rates depend on the electron density, the ionization balance of all elements 
are coupled. 

A rigorous treatment requires the electron density as an additional vari- 
able. However, the rate equation is then no longer linear and the rank of 
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the system becomes so large, that it is not feasible anymore for numerical 
computation. 

Furthermore, the radiative rates depend on the radiation field, which in 
turn depends on the occupation numbers. The radiation field dependence 
dominates the dependence on the electron density. Therefore, it is preferable 
to not treat electron density changes directly. In this work the consistency 
between the occupation numbers, the radiation field and electron density is 
obtained iteratively. 



2.8.3 Departure coefficients 

The departure coefficient b of a atomic level is defined as the ratio of the 
NLTE (or true) occupation number and the LTE occupation number 

bi = -g^ (2.73) 

for any level i in any ionization stage. In literature this is called the Zwaan 
definition [WZ72]. In the original Menzel definition [MC37], the departure 
coefficients 3 b* are not based upon the LTE occupation numbers, but on 
occupation numbers n*. These are computed per ionization stage with a 
LTE Saha-Boltzmann distribution based on the true NLTE continuum, i.e. 
equation (2.23) with n^J E replaced with n c . The subscript s for the ioniza- 
tion stage in equation (2.23) is omitted here. For the uppermost continuum, 
the bare nucleus, n* is defined as the true NLTE occupation number. 



n, = n c — exp 

g c 2 \2imi e kT 1 



Xi Xc 



(2.74) 



kT 

n* c = n c (2.75) 

with c being the continuum state to which the actual level i ionizes, and C 
the uppermost continuum. The b* are thus given by 

b* = ^ = b£z— = h/b c (2.76) 



nj n, 



c 



except for the last continuum, in which case the departure coefficient is 
converted as 



n. 



LTE 

C 



b* c = l = b c ^~ (2.77) 

n c 



3 here denoted with a * for notational convenience 
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2.8.4 Ionization balance 

When the rate matrix equation (2.72) is solved the ionization balance is 
determined. The ionization balance is important, as it determines the elec- 
tron donation of each species. In TE and LTE changes to the ionization 
balance exactly follow from changes to the (local) temperature via the Saha- 
Boltzmann law. In NLTE the influence of the temperature to the ionization 
balance is indirect, via the radiation field, and the rates. It cannot be de- 
termined directly. When changing the temperature, an assumption must be 
made about the change to the ionization balance, and therefore about the 
tii for all levels of all ionization stages of all elements. 

For this assumption two extreme approximations can be made. In the 
first approximation the n\ are completely independent of the local tempera- 
ture T . The rti are kept fixed, the LTE energy level distribution is computed 
for the new temperature, see page 2.3.2, and from that the bi follow. With 
the fixed ionization balance also the electron density becomes independent 
of the temperature. 

In the second approximation the rii adapt as much to the temperature 
as they would in LTE. The bi are kept fixed and the n; follow from the 
new LTE distribution for the new temperature. In this approximation it is 
convenient to directly compute the electron density from the bi, for which a 
NLTE partition function can be introduced, in analogy to the TE partition 
function (equation (2.13)). The total number density per ionization stage s 
is 

„LTE 

N * = E = E b ^ E = — ^ LTE ( T ' b *) ( 2 - 78 ) 

90s 

where i is the excitation state of ionization stage s and Os denotes its ground 
state. The NLTE partition function Q^ LTE (T) is defined as 

Q^ LTE (T, bi) = £ b ls9ls eM-Xis/kT) (2.79) 

i 

From equation (2.78) the number density of the ground state of ionization 
stage k can be written as 

n 0fe = J^gokb ok e- Xok/kT (2.80) 

Using equation (2.78) for s = k + 1 and the Saha-Boltzmann distribution 
(equation (2.17)) it follows 

V N k ) NLTE -Q^N c { V ) (A81) 

This equation is similar to the Saha equation (2.15) for TE. The ionization 
balance is solved analog to TE, which is described in section 2.3.2. 
The impact of these approximations is discussed in section 4.4 
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2.9 The spherically symmetric ID radiative trans- 
fer equation 

The time-independent, Lagrangian radiative transfer equation (2.4) in the 
case of a spherically symmetric atmosphere (SSRTE) can be written as 
[MW84] [HB99] 

ar ~dr + ° M ~dji + " A dA + 4axIx = r,x ~ XxIx ( 2 - 82 ) 
in which 

a r = 7 (^ + /3) 
atp = 7(1 - 



1 + 2/ , R sW 

7 + 



"a = 7 



with velocity v , /3 = v/c and 7 2 = 1/(1 — ^ 2 ). 

On the right hand side of equation (2.82), r]\ = rj\(r) and xx = Xx{f) 
are given by equations (2.24) and (2.25). 

The basic problem of radiative transfer is that the evaluation of a par- 
ticular I\(t\,(i) requires J\, and therefore I\ in all directions, and therefore 
the source function S at many locations and many wavelengths, and there- 
fore J on these locations and wavelengths. With equation (2.24) the SSRTE 
becomes an integro-differential-equation. In this work a monotonic velocity 
field is assumed. Then the SSRTE becomes a boundary value problem in 
spatial coordinates and an initial value problem in wavelength. The SSRT 
equation (2.82) is solved with the operator splitting method as described in 
[Hau92] and [HB04]. 

2.9.1 The Operator splitting method 

The mean intensity J\ is obtained from the source function S by a formal 
solution of the SSRTE, which is symbolically written using the A-operator 
A\ as 

Jx = Aa^a (2.83) 

In the following the subscript A is dropped to ease the notation. The usual 
Lambda iteration method 

Jncw = AS'old (2.84) 
Snew = (1 - e) J„ ew + e£ (2.85) 
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fails in the case of large optical depths and small e (thermal coupling pa- 
rameter, B is the Planck function). The range of the Lambda operator is 
only in the order of At ~ 1. At large optical depths, the mean intensity 
calculated with the Lambda iteration must be J\ = B\ + 0(e~ Tx ) and the 
convergence is very poor. A physical description of this effect can be found 
in [Mih80]. 

Therefore, the operator splitting method with an approximate lambda 
operator A* (ALO) is used, which is defined by 

A = A* + (A - A*) (2.86) 
and equation (2.83) is rewritten as 

J new = A*S ncw + (A - A*)S old (2.87) 
Using equation (2.85) it follows 

[1 - A*(l - e)] J„ew = Jfs ~ A*(l - e) Joid (2.88) 

where Jf s = AS'oid is the formal solution. With equations (2.86) - (2.88) new 
values of the mean intensity J ncw can be obtained and with equation (2.85) 
the new source function can be calculated. 



2.9.2 The Approximate Lambda Operator 

The formal solution is performed along characteristic rays [OK87]. Along 
the characteristic rays the SSRT equation (2.82) has the form [Mih80] 

Q- + a x — = r]-(x + 4a x )I (2.89) 

where s is the path length along a ray. Before the SSRTE can be solved the 
characteristic rays have to be calculated. The source function is interpolated 
piecewise linearly or parabolically along each ray. For the specific intensity 
I(Tj) the following expressions are obtained along a ray k 

I k (T t k ) = I k (rt 1 )eMT?-i-T l k )+ I S(r)eMr-r k )dr (2.90) 



I^ll.expi-Ar^ + Alf (2.91) 

T k is the optical depth along the ray k with i as the running index of the 
optical depth points. S is a generalized source function which contains 
terms of the wavelength derivative. It is defined in [Hau92]. With T\ = 
and T k _ x < T k it follows for the calculation of T k : 



Ar£i = (*i-i + *i)|sti-*il/2 



(2.92) 



2.10. RADIATIVE EQUILIBRIUM 



25 



Xi = Xi + 4aA,i is the effective extinction coefficient at point i and \s^_i — 
is the path length between point i and i — 1 along the ray fe. 
For A/k the following expression applies 

AI? = aj5i_i + + (2.93) 

a^, /3^, and 7^ are interpolation coefficients. The expressions for the coeffi- 
cients can be found in equations (23) to (25) in [Hau92] and in [HB04]. If S 
is interpolated linearly then is zero. For parabolic interpolation all three 
coefficients are non-zero. 

The choice of the approximate A-operator is very important to achieve 
rapid convergence. The ALO used in PHOENIX, a generalization of the tri- 
diagonal operator of [OK87], is given in [HSB94]. It is in some sense exact, as 
its elements are exactly the elements of the A operator. The approximation 
that is used is that only band diagonal elements of the full A matrix are 
used in the construction of A*. The expressions for the elements of A are 
given in equations (27) to (32) of [Hau92]. 

The fastest convergence of equation (2.88) is achieved with A* = A, but 
to construct the full A is very time consuming. In [HSB94] the operator 
splitting method with band-diagonal ALO operators is described and an- 
alyzed. In general, a ALO with a larger bandwidth will lead to a higher 
convergence rate, and therefore less iterations are needed to reach a given 
accuracy. But the ALO construction takes more time as more coefficients 
need to be computed. The optimal value for the bandwidth of A* depends on 
several conditions, described in [HSB94], like the computer architecture used 
and the type of model atmosphere that is computed. However, the strongest 
dependence is the quality of the initial guess to the solution. Therefore, the 
optimal ALO bandwidth is to be determined experimentally. 



2.10 Radiative equilibrium 

If an atmosphere is in radiative equilibrium (RE) then there are no net 
sources or sinks of radiation in the atmosphere and the radiative flux through 
each shell of the atmosphere is constant. Usually a target luminosity L is 
specified or an effective temperature T e g and a reference radius i?* . The RE 
condition in the Euler frame is 

4vrr 2 F obs (r) = 4vrr 2 J tff s (r)dA 

= 47rr 2 F obs (r) = 47ri? 2 J B x (T cS ) dX ^ 
= L = 47ri2 2 <7sBT e 4 fr 
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where H ohs is the Eulerian flux, i?Q bs is the target flux, related to the Lu- 
minosity L. In last step the Stefan- Boltzmann law has been used, and 

9-7T 5 t 4 

(2 ' 95) 

is the Stefan-Boltzmann constant. 

The RE condition in the Euler frame may also be written in differential 
form 

d(r 2 H ohs ) 

'- = (2.96) 

dr 

In the Lagrange frame this becomes [HBBA03] [MW84] 

^ + + —r 2 (J - K) + A\J + K + 2PH) = (2.97) 

Or Or r Or 

Using the first angular moment of the radiation transport equation (2.4) one 
obtains from both (2.96) and (2.97) an expression that is valid in both the 
Euler and the Lagrange frame 



/ 



fa - xx-h) dA = (2.98) 

In the Eulerian case rj, x and J are replaced by the Eulerian 7? obs , x obs and 
J obs . 

2.10.1 Model convergence criteria 

The first step in the modeling process is to make an initial guess for the 
atmosphere structure {T, n^} as function of the radius r. RE is not imposed 
at this stage. Subsequently, the atmosphere is corrected for deviations from 
RE iteratively using a temperature correction method (see section 2.11). 

Inserting rjx and x\ from equations (2.24) and (2.25) into equation (2.98) 
the scattering term a cancels out 



/ 



(fa ~ °xJx) ~ (Xx ~ °x)Jx) dA = (2.99) 

so that the relative deviation from RE can be written as 

_ f(r]\ ~ X\J\) dA 

£RE,dif = -ft \ja\ (2.100) 

J (xx — &x)Jx dA 

This is the differential form of the relative RE error. The definition is made 
in the Lagrange frame, because that is the frame the opacities are computed 
in. 
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The integral form follows from equation (2.94) 




— £RE,int 



(2.101) 



Note that this definition is made in the Euler frame. Defining it in the 
Lagrange frame is almost equivalent. But the target flux is formally specified 
in the Euler frame and the Lagrangian Ho is derived from it. Therefore, the 
H$ hs is the more fundamental quantity to check RE against. 

If the deviations from RE get below some threshold, typically £re < 1%, 
then the model is called converged. However, this condition is not very strict, 
since errors in single layers, as well as errors in the innermost or outermost 
regions (optically thick and thin respectively) are not as important as broad 
range errors in the middle region, that has the most important influence on 
the resulting spectrum. 



2.11 The temperature correction method 

One way to obtain radiative equilibrium is to use equation (2.100) in a sim- 
ple linear Newton-Raphson procedure. In the first step the temperature is 
changed arbitrarily. Subsequent temperature corrections can then be com- 
puted according to the induced changes in £re- However, this scheme is not 
very robust. It takes a lot of steps to converge to the required precision or 
it does not converge at all. 

A different approach is the Unsold-Lucy (UL) temperature correction 
scheme [Uns55] [Luc64]. In this method a direct correction AT to the tem- 
perature T is computed. This method is found to be much more stable than 
the Newton-Raphson scheme [HBBA03]. 

In this section all unprimed quantities are Eulerian and the primed quan- 
tities are Lagrangian. 

2.11.1 The Unsold-Lucy method 

The original procedure was developed for plane-parallel atmospheres. The 
derivation of the spherically symmetric generalization begins with the time 
independent, static, spherically symmetric transfer equation with an approx- 
imate LTE source function [HBBA03]. From equation (2.82) with /3 = 
follows 




l-p 2 d 
r dfi 



h = V\~ Xxh 



(2.102) 
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The zeroth and first angular moments are 

= >- 2 xa(Sa - Ja) (2.103) 

_^ = ^ + ?«iZ* ( 2 ,„4) 
Or r 

Integration over all wavelengths and assuming the approximate expression 
for the source function (see section 2.4) 

S x = (k x B x + a x J x )/(k x + a x ) (2.105) 

yields 

dr 2 H 

= r 2 ( KB B - kjJ) (2.106) 
dK 3K — J 

- XH H = — + 2.107 

Or r 

Here B, J and H are the integrals in wavelength of B x , J x and H x and 

KB = B I KxBxdX (2 - 108) 



= jj Kx J x d\ (2.109) 
XH = ^JxxH x d\ (2.110) 

are the mean opacities. Using the Eddington factor 

f = K/J (2.111) 

and multiplying equation (2.107) with the sphericity factor, introduced by 
[Aue71] 



r 2 f r Sf-l , 
q = -2 exp / - dr 



and writing ^ = r 2 i?, J? = r 2 J and Jf 7 = r 2 H one obtains from equations 
(2.106) and (2.107) 

KB^-KjJ? (2.112) 



dr 
dqff 
dr 



= -qXH-^ (2.113) 



It is assumed that the opacity means kb, kj and xh, as well as the Edding- 
ton factor / are independent of the temperature. Accordingly they do not 
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change from iteration to iteration. Equations (2.112) and (2.113) for the 
current iteration can then be subtracted from those for the next iteration: 

= k b A.% - kjAJ (2.114) 

dqf A J 

qJ Q / = -qxnAJf (2.115) 

with AJt? = J^ n ew — J^oid and J^ncw being the target flux and Jt? id the flux 
in the current iteration. Equation (2.114) is solved for ASS and equation 
(2.115) is integrated from the outermost layer R to r 

ASS = — ( _ + —A ^ (2.116) 

kb \ or or J kb 



rr ^ 

+ / q(r')xH(r')A^f(r')dr' 
Jr 



(2.117) 



IR 

The second Eddington factor 

g = J/H (2.118) 

can be used to express the A ^ ' (R) term in equation (2.117) in This 
factor is assumed to be independent of the temperature for the outermost 
model layer at radius R. Inserting equation (2.106) (for the J^ \d term) and 
(2.117) in (2.116) and dividing by r 2 yields 

AB(r) =AB 1 (r) + AB 2 {r) + } & f Q (2.119) 

r z Ks or 

ABi(r) =—J{r) - B(r) (2.120) 

K B 



AB 2 (r) 



k b r 2 q(r)f(r) 

(g(R)q(R)f{R) [J? (R) - JP(R)] (2.121) 
+ J\(r')XH(r') [Jfoir') - W(r')\ dr 1 

where the target flux is written as Jif new = Ji?o- 

A static atmosphere was assumed so that the radiative equilibrium con- 
dition implies a constant flux ^o(r) = J^o or equivalently (d,3%o)/(dr) = 0. 

All quantities at the right hand sides of equations (2.120) and (2.121) 
are available upon completion of each temperature correction iteration, so 
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that AB can be determined. The temperature correction AT follows from 
AB using the Stefan-Boltzmann law (equation (2.94)) 

AT(r) = -AB = - — (2.122) 

In the inner regions J\ ~ B\ so that the AB\ term becomes small. In a 
radially extended atmosphere, the 1/r 2 factor makes in the AB2 term small 
in the outer regions. 

2.11.2 The non-static generalization 

The non-static generalization of the UL temperature correction method must 
be derived from the complete SSRTE (2.82). 

However, if the expansion velocities are small < 1) then the result 
for the static atmosphere (equations (2.120) and (2.121)) can be used as a 
reasonable approximation when the target flux for the Lagrange frame H' is 
used. For an atmosphere in radiative equilibrium the flux in the Lagrangian 
frame is not constant, see equation (2.97). The value of the Eddington flux 
H , as measured by an external observer, is a given model constraint. The 
following transformation formulae relate the Lagrangian (primed quantities) 
and Eulerian systems for the wavelength integrated moments of the radiation 
field [MW84] 

J = 7 2 (J' + 20H' + (3 2 K') (2.123) 
H = 7 2 ((1 + p 2 )H' + (3( J' + K')) (2.124) 
K = 7 2 (K' + 2pH' + /3 2 J') (2.125) 

The Eulerian target flux Ho can be converted to the Lagrangian, using the 
inverse of equation (2.124). 

H' = 7 2 H((l + /3 2 )-f3(j + k)) (2.126) 

where j = J/H and k = Kj H. If j and k are independent of the temperature 
then j ncw = j id an d k ncw = A; id. For the next iteration (denoted with new) 
equation (2.126) can be written as 

^ = 7 2 #o((l + /3 2 )-/3^) (2.127) 

Here the target fluxes Hq and Hq are inserted H new = Hq and H' new = Hq. 
J, H and K are the values obtained in the current iteration (denoted with 
old). 

This equation yields the target Lagrange flux H'(r) for all radial points 
r in the atmosphere. The temperature correction for moving atmospheres 
is then given by equation (2.122) using the Lagrangian variables J', H' and 
K' . 



Chapter 3 



Nova atmosphere 
construction 



In the theoretical description of nova outbursts, described in section 1.1, 
multiple phases can be discerned. Whereas the agreement between theory 
and observations in the optical and UV was very good, the X-ray observa- 
tions do not fit into the picture. This requires further investigation in order 
to improve our physical understanding of the phenomenon "nova" . It was 
described in that same section that X-rays gradually start to emerge from 
the atmosphere after the "fireball" phase when the optical lightcurve has 
already been steeply declining for a while. The phase identified with strong 
X-ray luminosity is called the "supersoft source (SSS) phase". It is this 
phase that is modeled in this work. 

According to theory [Sta89] the density profile in earlier stages can be 
described by a power law 

p(r) oc r~ N (3.1) 

In early phases the parameter N is of the order of 15. In later phases the 
value of N becomes much smaller, around 3. Successful modeling of early 
optical and UV lightcurves were obtained with PHOENIX based on density 
structures of this type [HSS + 95]. 

In the SSS phase the total optical thickness of the envelope has decreased 
further, so that the radiation emerges from deeper regions in the extended 
nova atmosphere. It was described in section 1.1 that the part of material 
that is accelerated in the outburst below the escape velocity falls back onto 
the white dwarf atmosphere and reestablishes hydrostatic equilibrium. In 
the SSS phase a stellar wind-type mass loss occurs, driven by radiation 
pressure on atomic lines [Sta89] . This is observationally confirmed by finding 
systematic blueshifts in absorption line centers [NSB+03] and [NSB+07]. 

Line driven stellar winds are theoretically described by a beta law, see 
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[CAK75], [KPPA89], and [Lam97]. 

v(r) = v* + (voo - v*) (l- ^) (3.2) 

with being the velocity at the bottom of the wind r* , and the terminal 
velocity. The value for f3 obtained for hot giant winds from theory and 
observations is /3 = 0.8. For cool stellar winds the values determined from 
the observations are typically larger [Baa92], [BKR + 96]. 

It cannot be known a priori how optically thick the wind is that escapes 
from the hydrostatic white dwarf atmosphere. Therefore, in this work a 
hybrid-type atmospheric structure is assumed, composed of a hydrostatic 
base and an expanding envelope. This is based on an idea of [AufOO]. In 
the models the optical thickness of the wind can be varied freely and the 
hydrostatic core provides proper boundary conditions at the base of the 
wind. 

In the rest of this chapter the subsequent steps involved in the atmo- 
sphere construction procedure for the hybrid-type nova atmospheres are 
described. 

3.1 Atmosphere construction procedure 
3.1.1 Envelope: Density 

Following [HSS + 95] the envelope is considered expanding but stationary, 
and spherically symmetric. Therefore, it is assumed that all time dependent 
terms both in the hydrodynamics and in the radiative transfer equation can 
be neglected and all quantities depend only on the radial coordinate (except 
the specific intensity I\ of the radiation field which, in addition, depends 
on the angle to the radial direction and A). These assumptions are justified 
since the hydrodynamic timescales are much larger than the timescales for 
photon thermalization and for the establishment of the statistical equilib- 
rium (see section 2.8.1). The assumption of spherical symmetry is only an 
approximation since novae are contained in a binary system which is clearly 
not spherically symmetric. Whereas the 1-dimensional treatment presented 
in this work is already very computationally expensive with the currently 
available computer power, the requirements for a 3-dimensional generaliza- 
tion are hardly met by the largest supercomputers presently available. Fur- 
thermore, although 3D radiation transport is doable with PHOENIX [HB08], a 
3D 'stellar' structure is needed as input, and at present no self-consistent 3D 
hydrodynamical nova evolution models are known in literature that could 
provide such a structure. Until such input structures become available the 
most realistic ID models can provide the best possible constraints that are 
important for the development of nova evolution theory, which in turn is the 
basis for detailed 3D hydrodynamical models. 
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Since the wind is assumed to be time-independent and stationary, the 
rate of mass loss is constant with time and therefore with radius 

M(r) = 4vrr 2 p(r)v(r) = M (3.3) 

This is called the continuity equation and represents mass conservation in 
the envelope. With the velocity field being prescribed by equation (3.2) the 
density structure of the envelope is fixed 

4irr 2 v(r) 



3.1.2 Envelope: Radial sampling 

The radial sampling of the atmosphere should provide an optimally smooth 
grid for the source function to be interpolated over, see section 2.9.2. Large 
steps in optical depth from one layer to the next when integrating the spe- 
cific intensity along characteristic rays leads to inaccurate solutions for the 
radiation field, see section 4.5. Improper radial sampling can in principle be 
remedied by increasing the number of layers the atmosphere is divided in, 
but that is computationally expensive. So in order to optimize the compu- 
tation time and accuracy of the radiation transport the radial sampling is 
important to consider. A number of radial sampling prescriptions based on 
power laws or the cosh function can be found in literature (e.g. [SSMS97]). 
But the sampling rates obtained from these prescriptions are far from ideal. 

The radial extension of the envelope can be very large (typically a factor 
100-1000) and the opacity scales roughly with the density. Therefore, the 
radial grid is constructed thus that the density is divided into logarithmically 
equally spaced parts. 

p(l) = p*(Ap) 1 *- 1 (3.5) 

where I = 1,2,...,/* is the layer number increasing inwards, and Z* the 
bottom layer of the wind (and simultaneously the top layer of the hydrostatic 
core). The determination of Ap is possible because the density structure is 
fixed by the parameters of the wind, equation (3.4). Now the radial grid 
follows from the continuity equation (3.3) 

= «'» 2 " ») (3 ' 7) 

which is straight forward to solve numerically, layerwise, starting from the 
inner boundary, for example with a Newton-Raphson scheme. 



34 



CHAPTER 3. NOVA ATMOSPHERE CONSTRUCTION 



3.1.3 Envelope: Optical depth scale 

Once the density and radial grids for the expanding envelope are fixed, the 
optical depth grid 1 T\ icf for that range (/ = 1,...,/*) can be integrated 
inwards using 

dr A rof = - (XA ref - xa Xlci ) dr (3.8) 

where x is defined by equation (2.25) except that the cores of atomic lines 
are excluded. It is important to include the opacities of far line wings in the 
construction of the t\ ic{ scale, since these make up for a significant part of 
the total opacities in some parts of the atmosphere. This has a major impact 
on the t\ tc{ scale, x is a free parameter < x < 1 that is used to turn off the 
scattering contribution to the extinction in the optical depth reference scale, 
which mainly influences equation (3.8) of the outermost layers. The values 
of x and A rc f depend on the physical conditions of the atmosphere. For the 
models computed in this work the typical values are x = 1 and A re f = 40 A. 

3.1.4 Core: Optical depth scale 

For the hydrostatic layers (/ = /* + 1, ... , N) the structure is computed on 
a fixed logarithmically spaced optical depth grid 

log(r(0) = log(r(Z*)) + M 

starting from the value of the dynamic region's bottom layer r(Z*). Here 
r = T\ ref to ease the notation. In the following for all occurrences of r 
referring to the optical depth grid r = t\ ic1 holds. The lower boundary is 
defined for an optical depth t(N) = r max with a typical value of r max = 1000 
to ensure it is set well below the thermalization depth 2 . 

3.1.5 Both: Temperature structure 

The temperature structure T(l) is provided as input to the model. Devi- 
ations from that temperature structure that fulfills the condition of radia- 
tive equilibrium are iteratively corrected for, see sections (2.11.1), (4.1) and 
(4.8.1). 

The effective temperature specifies the Eulerian target flux Hq according 
to equation (2.94) (see also equations (2.127) and (2.121)) via the Stefan- 
Boltzmann law 

O) = 5jU B T e 4 ff (3.10) 

1 This r-scale plays no direct role for the expanding envelope, but it provides the bound- 
ary condition for the construction of the hydrostatic core of the model (see section 3.1.4). 

2 The thermalization depth is defined as the optical depth at which the radiation field 
"thermalizes" J\ — B\ as seen from outside. It is not easy to pin down the exact location 
at which this is satisfied, therefore mostly a thermalization proximity is meant. 



log(r(iV))-log(r(^)) 
N-L 



(3.9) 
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In this definition the 'effective radius' (the radius specified together with 
the effective temperature in the Stefan-Boltzmann law) is not related to the 
optical depth scale but fixed to the outer radius of the white dwarf, being 
the base of the wind. This choice has a large advantage in the convergence 
properties and the stability of the models. The optical depth scale changes 
during the iterative process, so that also a target flux that depends on it 
would fluctuate. 

Care must be taken with the physical interpretation of the "effective 
temperature", since in fact it only directly specifies the model's luminosity 
(and not any real temperature). 



3.1.6 Both: Structure completion 

Now the complete atmospheric structure can be computed for each layer. 

For the expanding envelope the gas pressure P g (l) is calculated from the 
ideal gas law using the given density and temperature structures 

P S = ^ (3.11, 

fj,m H 

with \x = fh/niH being the mean atomic weight per free particle (including 
free electrons) in atomic units, and P g = ^ Pi = being the gas pressure 
is the sum of all partial pressures, where i specifies all types of particles 
including free electrons. 

For the hydrostatic core this is done by numerically integrating the hy- 
drostatic equation 



dP 9 _ 9 



dr st d K s td 



(3.12) 



downwards, starting from the transition layer at which the pressure, den- 
sity and thus the mass extinction coefficient are already known. Finally, the 
radius for every layer is calculated from equation (3.8). 



3.2 Structure check 

Thus far no restrictions were made to ensure a smooth transition from the 
expanding to the hydrostatic region. The hydrostatics is parameterized by 
the gravity log(g). The wind is parameterized by mass loss rate M, the 
terminal velocity i>oo, the wind parameter j3 and the base velocity i>*. 

It is right this last parameter v* that is responsible for a gradual transi- 
tion in the boundary layer /*. The effect of this parameter on the transition 
is shown in figures 3.1 to 3.3. Apparently, a bad value for vq distorts the 
structure of the core. Whereas in a hydrostatic atmosphere the density 
(layer-)gradient (figure 3.2) is constant, this is not the case if the boundary 
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Figure 3.1: The optical depth scale, here plotted against the density, shows the transition 
from core to envelope around p — 10~ 5 g/cm 4 as a kink in the curve. The bottom wind 
layer is marked with a cross sign, the other layers with a plus. 

The base velocity of the expanding envelope v* has an important impact on the transition. 
The proper value (red curve) yields a curve that is straight at both sides of the intersection 
point (transition layer). If the wind starts too far 'in' or 'out' (black and green) then 
the boundary condition for the hydrostatics is inappropriate, and the hydrostatics are 
distorted. 



conditions supplied by the wind, i.e. the opacity and density of the transi- 
tion layer, are not appropriate. This distortion can have significant influence 
on the structure and thus also on the resulting spectrum. Another effect of a 
unsuitable value for vq is that the atmosphere cannot be accurately brought 
into radiative equilibrium. 

As a last step in the atmosphere construction procedure, this transition 
is checked, and corrected if too far off. This is done by checking if the density 
sampling rate (see figure 3.3) is approximately constant for the hydrostatic 
region. If dp/dl at the outermost hydrostatic layer is off from the mean over 
the 10 inner layers (excluding the innermost, with optical depth of r > 10 2 , 
as these are not reliable) by more than 20%, then the value of is corrected 
and the atmosphere construction scheme is repeated. 
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Figure 3.2: Here dr le t/dp, being the slope of the graphs in figure 3.f, is plotted log- 
arithmically against the layer number. The discontinuity at layer 67 is inherent to the 
construction with an ad-hoc boundary between hydrostatics and hydrodynamical expan- 
sion. The bottom wind layer provides the boundary condition for the hydrostatic core. 
The parameter «* must be tuned in order to provide the proper conditions at the boundary 
that meet the requirements of hydrostatics. In this example «*/«oo = 2 • 10 -4 (red curve) 
is the proper value. With too low or too high values (black and green) the hydrostatics 
(on the left side of the transition) are distorted. Close to the transition, the lines get 
curved away from the straight, undistorted run for hydrostatics. 

In the undistorted case the slope of the hydrostatic curve only depends on the radial sam- 
pling rate (number of layers per optical length) and \og(g). In this example the \og(g) is 
equal for the three models and the different slopes are caused by different sampling rates 
due to a shift 'out'-, or 'in'-wards of the transition layer corresponding to the value of v*. 
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Figure 3.3: This plot shows the sampling rate of the density by the radial grid dp/dl. 
The sampling of the expanding envelope (I < 67) is fixed by construction (equation (3.5)). 
The sampling of the core (I > 67) follows from the hydrostatics on a logarithmic r re f scale 
with the transition layer as boundary condition. 

With the proper boundary condition (red curve), the density sampling rate for the hy- 
drostatic range is approximately constant. For too small values of v,/voo the density 
sampling dp/dl can become smaller than 1, in which case the density structure is no 
longer monotonically increasing. For too large values, the changes in density from layer 
to layer become large which is disadvantageous for the computation of the radiation field 
(especially, the interpolation of the source function over large optical lengths can introduce 
large inaccuracies, see sections 2.9.2 and 4.5). 



Chapter 4 

Methods for good NLTE 
models 

4.1 Global iteration scheme 

It is a complex task to attain a fully converged NLTE atmosphere model. 
The procedure used in PHOENIX is described in [HB99]. Basically, it is a four- 
steps iterative process. It starts from an initial guess for the atmospheric 
structure, something like [T, p,p gas , nj](r). 

1 In the first step the radiation field is computed using the SSRTE (2.82) . 

2 Then, from the radiation field the deviation from radiative equilibrium 
(equation (2.100)) can be checked and, if necessary, the temperature 
correction (equation (2.94)) can be computed. 

2b For NLTE models the radiative rates are calculated from the radia- 
tion field, and the rate matrix equation is solved. This updates the 
occupation numbers n{. From these Hi the bi are determined. 

3 The new temperature T ncw (r) is applied. 

3b For NLTE models the n« for this temperature follow from n^ TE (T ne w) 
and the bi (computed in step 2b), using equation (2.73). 

4 The structure is adapted to the new temperature. 

With the updated structure, after step 4, a new cycle starts. This is con- 
tinued, until in step 2 the deviation from RE is smaller than a prescribed 
threshold. 

In this work, the conventional iteration scheme used in PHOENIX was 
found to have a poor convergence performance for NLTE models. Often the 
solution diverges or very many cycles are needed to attain the prescribed 
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convergence criteria. This is sometimes referred to as apparent convergence 
as opposed to true convergence [PHNS05]. 

Relaxing the assumption of LTE to Non-LTE gravely complicates the 
coherence in atmospheric models. In LTE the atomic level populations of the 
gas are fully determined by one single quantity, the temperature T (section 
2.4). In such models the temperature structure is one of the most important 
properties of the atmosphere, it characterizes the spectrum. 

In NLTE the level populations are no longer directly determined by the 
temperature, but by the rates. In the outer parts of the atmosphere, where 
the densities are low, the radiative rates play a dominant role compared to 
the collisional rates. The bound-bound radiative rates are functions of the 
mean radiation field J only (equations (2.33) and (2.45)). The bound-free 
radiative rates primary depend on J, but also on the electron density n e 
and the temperature T (equations (2.65) and (2.66)). 

In order to improve the convergence performance for NLTE models the 
iteration procedure must be altered. It is important that the radiation field 
is fully consistent with the population numbers, before the temperature 
corrections are derived and the structure is updated. In the following this 
is called radiation-matter consistency. In LTE models, this consistency is 
immediately given. In NLTE models, for this consistency an inner iteration 
loop is needed, comprising of step 1 and step 2(b) only. Or, to put it 
another way, in every iteration the structure must not be updated (step 
3 and 4) before this consistency is given. In the following, the iterations 
in which the structure is not updated are called idle cycles. Whether this 
condition of (approximate) consistency between the radiation field and the 
level populations is given can be checked by the changes in the J and the 
rii from iteration to iteration. 

An implicit condition for this consistency is that the ionization balances 
of all NLTE species do not change significantly from iteration to iteration. 
Here, ions with small ion number densities are less significant than those 
with large. In order to check for this condition the following definition of a 
weighted average of ionization change is useful 



where the sums extend over all ionization stages s of chemical element E, N s 
is given by equation (2.78), and the max() function returns the maximum 
value of its arguments. 

In figure 4.1 an example is given of a model spectrum resulting from 
poor radiation-matter consistency. The model performed the same number 
of updates to the structure as the reference model, performed with idle 



AN E = 



£ N s max(N 8 /N 8j0ld , N 8j0ld /N 8 ) 
£max(iV s 2 /iV Si0ld ,iV s 

,old) 



(4.1) 
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Figure 4.1: In NLTE models the radiation field and the population numbers must be 
consistent before the atmospheric structure is updated. This consistency is achieved using 
idle cycles, in which the radiation field, the radiative rates and the population numbers 
are updated, but the structure is held fixed. Without idle cycles NLTE models hardly 
converge to RE or often even diverge. 

This plot shows two spectra, for which the same number of updates to the structure were 
performed. The models are identical except for using idle cycles or not. The model with 
idle cycles (red) is in RE, whereas the model without idle cycles (black) diverges (becomes 
worse after more iterations). 

cycles. The model without idle cycles diverges and the resulting spectrum 
becomes worse in every step. In the reference model the number of idle 
cycles per full cycle is one. 

Typically, one or two idle cycles are needed, depending on the changes of 
the radiation field. If the changes are small, the radiation-matter consistency 
is quickly attained. 

4.2 Pure NLTE 

NLTE models are commonly treated as a gradual refinement to LTE models. 
That means that an LTE model is used as the initial guess, on which more 
and more species are subsequently treated in NLTE. 

This approach may be good for models that are close to LTE in all re- 
gions of the atmosphere. But if the deviations from LTE are large, then 
LTE level populations are far off from the actual populations. Consequently 
LTE opacities are far off from the actual opacities. The opacities dictate the 
radiation field and the atmospheric structure. It is a matter of experimen- 
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stage 


atomic species 


# lines 


# levels 


1 


H He C N 


2800 


600 


2 


H He C N Ne Mg Al Si S Ar Ca 


46000 


4000 


3 


H He C N Ne Mg Al Si S Ar Ca Fe 


160000 


7600 



Table 4.1: Every model is computed in three subsequent stages with an increasing 
number of chemical elements treated in NLTE. The table shows the typical number of 
atomic levels and lines considered in each stage. 

tation to find out if, for the model under consideration, the LTE opacity is 
more realistic than a null-opacity (no opacity at all) for the specific species. 
This depends on the actual deviation from LTE. 

For the models performed in this work, it is found that adding LTE 
opacities is a worse approximation than leaving those specific opacities out 
altogether. This effect is shown by an example of a typical model in figures 

4.2 and 4.3. A comparison between pure NLTE and pure LTE models is 
made in section 5.2. 

The comparison is made between a full stage 2 model, which treats all 
elements up to Calcium in NLTE (see Table 4.1), and a model in which 
oxygen is not treated at all or treated in LTE. The two test models started 
from the reference model, and passed the same number of iterations. In all 
models RE is obtained. Clearly, the model without oxygen opacities is more 
similar to the full model than the one with LTE oxygen opacities. Note that 
in the test model with oxygen in LTE still the far majority of the opacities 
was treated in NLTE. 

Generally, adding LTE opacities to the type of NLTE models considered 
in this work detoriates the resulting structure and the model spectrum. 
Therefore, only pure NLTE models are computed, unless stated otherwise. 

4.3 Modeling steps 
4.3.1 Multi-stage approach 

The number of iterations needed to let a model converge to RE depends 
significantly on how close the initial guess is to the final model structure. 
The time needed for an iteration depends on which chemical species are 
treated in NLTE, especially the number of atomic lines (step 1) and levels 
(step 2b). In order to minimize the total computation time the modeling 
process is divided into several stages. In every stage the number of species is 
increased, see Table 4.1, and thus the model comes closer to the final result. 
This is a typical approach, see [PHNS05]. 

The first stage quickly delivers a rather good initial guess for the second 
stage. And the second delivers relatively fast a good initial guess for the 
last stage, being the final model, which is very computationally expensive. 
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Figure 4.2: Two test models (black) are compared to one reference model (red). In The 
reference model the following elements are treated in NLTE: H, He, CNO, Ne, Mg, Al, Si, 
S, At and Ca. It has no LTE opacities. In the test model in the upper graph oxygen is 
treated in LTE. In the lower graph oxygen is not considered at all, so no oxygen opacities 
are included. 

Clearly, the oxygen LTE model departs more from the reference model than the model 
with no O at all. Note that still the far majority of the opacities is in NLTE. 
In general, adding LTE opacities to NLTE models of the kind presented in this work 
deteriorates the results. 
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Figure 4.3: Two test models (red and green) are compared to one reference model 
(black). The models are described in the previous figure caption. The upper graph shows 
the temperature structure against the density. The lower graph shows the optical depth 
at a reference wavelength (40A) r rc f against the density. This r rc f is used to construct the 
atmosphere. 

Clearly, the oxygen LTE model departs more from the reference model than the model 
with no O at all. Note that still the far majority of the opacities is in NLTE. 
In general, adding LTE opacities to NLTE models of the kind presented in this work 
deteriorates the results. 
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An example of the resulting spectra and temperature structures for the tree 
stages are shown in figure 4.4. Note that in every step the model comes 
closer to the final result. 

4.3.2 Adding new species in NLTE 

If new species are added in NLTE, then the occupation numbers rij = n^ TE 
for those species are far off. That means that initially the opacities, therefore 
the radiation field, and therefore the radiative rates are far off from their 
right values. The rij of the NLTE species of the previous stage are already 
close to the converged solution. Keeping those rii fixed for a few idle cy- 
cles when adding new species in NLTE prevents them from being detoriated 
by the yet unrealistic new opacities. Furthermore, this improves the radi- 
ation field which helps the new species to converge to the final occupation 
numbers. 

When the n; of the new species have converged to an acceptable accuracy, 
i.e. when their changes from one idle iteration to the next become small, 
then the rii of the old species need to adapt to the altered opacities. The 
same trick is applied again, keeping the new rii fixed and let the old converge. 
After that, normal model iteration can proceed, as described in section 4.1. 

4.3.3 Countering micro-inconsistencies 

In NLTE models there is a huge number of variables, for each radial point in 
the atmosphere, e.g. temperature, electron density, and occupation numbers 
of all atomic levels. Using the iterative modeling process, see section 4.1, 
their coupling is treated in an indirect way, via the radiation field. 

Even with the greatest caution, it can occur that during the iterative 
refinement procedure some occupation numbers enter a wrongly or weakly 
coupled region. In such a case, once a (couple of) rii is in such a region it will 
not be able to approach the right value soon. In this work this is called micro- 
inconsistency, as generally their influence is small on wavelength integrated 
quantities, and thus on the atmospheric structure. Especially, if the initial 
structure is far from the final solution there is a long way for such micro- 
inconsistencies to develop. 

A rigorous way to deal with micro-inconsistencies is to sequentially reset 
parts of the variables to default values and let them iteratively converge 
again, keeping all others fixed. In principle, this can be done with any part 
of the variables, like the outer part of the temperature structure. But the 
temperature is directly coupled to the gas pressure and density via the ideal 
gas law. Therefore, it is easier restrict to a ni-reset method, resetting the rii 
of one or a few species at a time. As default values the LTE values can be 
used. A nj-reset is in fact a set of idle-cycles in which not only the structure 
is kept fixed, but also the n« of all but a small, subsequently changing set of 
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1Q -2 10 -4 1Q -6 1Q -8 1Q -I0 1Q -12 1Q -14 
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Figure 4.4: The modeling process is done in three subsequent steps (black, red and 
green), in which the number of atomic species is increased, and with that the computation 
time. In every step the model comes closer to the final model (green). This procedure 
provides a computation time saving way to obtain a converged complex NLTE model. 
The upper graph shows the spectra, the lower graph shows the temperature structures for 
the three subsequent stages. 
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Figure 4.5: In the process of iterative refinement of an initial structure to the final 
solution for the model micro-inconsistencies can occur. That means that some of the large 
number of rn in the model end up with a wrong value and diverge or slowly converge. 
Most often, micro-inconsistencies are not very noticeable, but can yet have a significant 
influence on the spectra. 

In the example shown here, the existence is noticeable: a strong emission line appears due 
to a wrong m (black curve). After a n^-reset (see section 4.3.3) the problem is solved (red 
curve) . 

the atomic species. The number of idle cycles needed per set of species to 
let the rtj to approach the right values depends on the model atom and on 
the deviation from LTE. 

Micro inconsistencies sometimes show up as strong emission lines, as in 
the example shown in figure 4.5. In such cases, their existence is evident. 
However, mostly they are inconspicuous but can still have a large overall 
effect on the spectrum. In any case, their existence can be checked for and 
eliminated by the n^-reset method. Therefore, in each of the three stages, 
a seemingly converged model is checked with the n,-reset method. After 
the reset, the model is iterated further until RE is reached. This process 
is repeated until the changes to the structure, and thus to the spectrum, 
between two n^-resets are small. 

4.4 Transferring temperature corrections to the n{ 

In the previous sections of this chapter 4.1, 4.2 and 4.3 it has been de- 
scribed, that bringing a model into statistical, radiative and hydrostatic (or 
dynamic) equilibrium is a numerically complex task, due to the large number 
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of indirectly coupled variables. In order to attain radiative equilibrium the 
temperature structure is subsequently improved, and after each correction, 
the rii must find their balance again. 

As described in section 2.8.4 the change in the local NLTE ionization 
balance due to a change of the local temperature can not generally be pre- 
dicted. In the limit of small departures from LTE, the NLTE ionization bal- 
ance changes according to the Saha-equation. In the opposite limit of large 
deviations from LTE, no coupling to the local temperature is expected. 

Ordinarily, in PHOENIX, when the temperature structure is corrected, the 
b* are kept fixed, and the new rii follow from equations (2.74) to (2.76) after 
solving the NLTE ionization balance, equation (2.81), for the new tempera- 
ture. This is a good approximation for the inner regions of the atmosphere. 
But analysis shows that in the outer regions this often leads to overcorrec- 
tions, from which the rii partly recover in following idle iterations, and finally 
approximate the balance before the overcorrection was made. An example 
of this effect is given in figure 4.6. These overcorrections needlessly desta- 
bilize the rii and therefore increase the number of idle iterations needed to 
stabilize it again. Also, this increases the problem of micro-inconsistencies. 

Also, the rii overcorrection affects the temperature correction method. 
The method is based on the assumption, that the ratios of the integrated 
opacities are insensitive to the temperature, equations (2.112) to (2.115). In 
the outer regions, where the densities are low, the contributions of strong 
lines become more significant in the opacity integrals. If the rii of these 
lines are destabilized by temperature corrections, then the opacity integrals 
are found to become sensitive to the temperature, and the method becomes 
inaccurate. In figure 4.7 the resulting temperature structures are shown for 
a model in which the rii are scaled with temperature corrections the ordinary 
(unattenuated) way, compared to a more suitable attenuated way, which is 
described below. 

In order to improve the way the temperature corrections are transferred 
to the rii the "small departure from LTE" approximation must be alleviated 
for layers in which the departures are large. For that purpose an approximate 
measure is needed for the departure from NLTE for each chemical element 
E that is treated in NLTE. One obvious way to define this departure is to 
average the bi of each element, weighting each state i by its number density. 

r EMi Y.K^rii 
oe = = — ^ 

But since the bi can be very large or very small numbers they dominate 
the sum, by the means that unimportant states can dominate the average. 
Therefore, it is better to average the logarithms of the bi 

OE = ^ = ^ 

E rii E rii 



(4.2) 



(4.3) 
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128 96 64 32 

Layer 

Figure 4.6: Ordinarily, the population numbers are scaled in a LTE-like manner with 
the temperature, keeping the bi or the b* fixed. That often leads to strong overcorrections 
in regions where the m are weakly coupled to the local temperature. 

The upper graph shows the ionization balance of oxygen, with the deepest atmosphere 
layers at the left. The solid curves show the NLTE and the dashed curves the LTE balance. 
In the middle graph, the relative changes from iteration 1 to iteration 2 are plotted on a 
logarithmic scale. The changes in the LTE balance show that the temperature changed. 
The changes in the NLTE balance follow the LTE changes, but are not identical since not 
the bi but the 6* were kept fixed. Obviously, overcorrections from iteration 2 in the outer 
40 layers are largely canceled again in the third iteration (bottom). 
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Figure 4.7: If the NLTE occupation numbers in the outer regions of the atmosphere 
are assumed to depend on the local temperature like the LTE numbers, then a number of 
problems are introduced. One of the problems is that the temperature correction method 
becomes inaccurate, which leads to wrong temperature structures. 

This plot shows the temperature structures of two models in comparison with one original 
structure (black) that was used as starting point for the models. The model with the 
ordinary full LTE correction to the m is shown in red. In the other model (green) the 
corrections were attenuated by a factor of 10 (xe = 0.1) in regions where the fa; were large, 
as described in the text. 



This measure 1>e is then used to attenuate the full, LTE like, corrections to 
the Hi. The rij can then be computed using a linear combination of the two 
limiting cases 

n 4 = x E [b in l TE (T new )] + (1 - x E ) [hnl TE (T old )] (4.4) 

where xe is the attenuation factor for chemical element E. There are unlim- 
ited possibilities to relate 1>e to xe, with the restriction that in the limiting 
case of 1>e ~ (note the logarithm in the definition) the attenuation must 
be off: xe = 0. 

While experimenting with relations between i>E and xe it was found that 
using the very simple xe = globally, i.e. the limit of no thermal coupling of 
the population numbers for all regions of the atmosphere, gives good global 
convergence. In those regions of the atmosphere where the departures from 
LTE are small the are found to quickly converge to their new balance 
for the new temperature. In contrast, in the regions with large departures 
the rii converge slowly. Using the approximation that is well suited for the 
complex regions still leaves enough time for the simpler LTE like regions to 
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converge, even if the n^-adaption approximation is bad. The performance 
of this approximation is shown in figure 4.8, which is the same model as in 
figure 4.6. In figure 4.9 the weighted averaged density corrections ANe, as 
defined by equation (4.1), are compared between the fully thermal coupling 
approximation and the no thermal coupling approximation. In the inner 
regions the thermal approximation is a bit better, but in the outer regions 
the non thermal approximation is far superior. 



4.5 Interpolation of the Source function 

In order to compute the specific intensity along a characteristic ray, see 
section 2.9.2, the source function is integrated, equation (2.90). The source 
function is only computed on the intersection points of the ray with the 
concentric shells, here often called layers, of the model. For the integral 
an approximation of the source function between the intersection points is 
required. If S is interpolated linearly or parabolically, then the integral can 
be written as a linear combination of S on either two or three intersection 
points, equation (2.93). 

It is not possible to generally favor the linear or the parabolic interpo- 
lation based on physical arguments. In PHOENIX a threshold optical depth 
value can be set. Below the threshold, for small optical depths along a char- 
acteristic ray, the linear interpolation is used, and above the threshold the 
parabolic. This allows to compare results with the linear and parabolic in- 
terpolation methods. The reason for the threshold switch method is that the 
parabolic interpolation occasionally gives very unrealistic values, so that the 
formal solution for the characteristic ray becomes inaccurate. This leads to 
inaccurate mean intensities which not only degrade the convergence proper- 
ties of the operator splitting method, but also lead to an inaccurate solution 
of the whole radiation field. These unrealistic interpolation values usually 
occur in small optical depths along rays. Setting the threshold high enough 
enables to avoid these problems. 

Comparative tests performed in this work showed that the linear and 
parabolic S'-interpolation method can lead to very different solutions of the 
radiation field. In the search for a more sophisticated method, that should 
have been able to determine the preferred interpolation method for each 
individual situation, a different approach was found. 

Plotting the source function at the intersection points against the op- 
tical depth along a ray, shows a contrarily tendency of both interpolation 
methods: if the linear method overestimates S between two points then 
the parabolic method underestimates it, or opposite. An example of these 
tendencies is given in figure 4.10. The plot shows the intersection points 
as black crosses, the linear and parabolic interpolations and a combined 
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Figure 4.8: In this plot the same model is shown as in figure 4.6, but now the non- 
thermal limit approximation was used for the whole atmosphere. From the first to the 
second iteration the temperature changed, but no direct changes are induced in the NLTE 
ionization balance. In the third iteration the changes to the balance in the outer regions 
(layer numbers smaller than 40) are still small. The strong overcorrections in these layers, 
as seen in figure 4.6, are avoided. This not only saves idle cycles, but is also very important 
for attaining a consistent set of population numbers. 
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Layer Layer 

Figure 4.9: The weighted averaged density corrections ANe, as defined by equation 
(4.1), are compared for oxygen in five subsequent iterations between the fully thermal 
m coupling approximation (left) and the non thermal approximation (right). In the first 
iteration the temperature was corrected, and in the fully thermal approximation the m 
were corrected accordingly. In the non thermal approximation the rii are not corrected for 
the temperature change. The changes visible in the top right plot are due to gas density 
changes in the atmosphere setup. 

In the fully thermal approximation even in the third iteration after the temperature cor- 
rection the ANe for oxygen is still 30% in the outermost layers. In the non thermal 
approximation the rij of the outermost layers are corrected only by 10% in total. In the 
lower layers the fully thermal approximation is superior, but the corrections are small and 
not problematic in either case. 
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linear 
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combined 




Figure 4.10: The source function along a characteristic ray is only computed on intersec- 
tion points of the ray with the concentric shells of the model. When the source function is 
integrated along the ray, in order to compute the formal solution of the radiative transport 
equation, the run of the source function between the points is needed. 
The plot shows an example of the source function against the optical depth along the ray 
on a double logarithmic scale. The intersection points are marked with black crosses. If 
the linear interpolation (blue) overestimates the source function then the parabolic inter- 
polation (red) underestimates it, and vice versa. A linear combination of the linear and 
parabolic case, here to equal parts (equation (4.5)), combines these tendencies and gives 
a much better approximation (green). 



interpolation. The combined case shows 

S(r) = (S hn (r) + S P Ur))/2 (4.5) 

Since the integral in equation (2.90) can be written as linear combina- 
tion of S on a few intersection points, it can also be written as a linear 
combination of the solutions with both interpolation methods. 

AI k = (gojjjg + fe <ar)^-l + «4in + ^par)^ + &7*pAl ^ g) 

a + b 

This method has the advantage, that the over- and underestimating tenden- 
cies of the linear and parabolic method are canceled for good parts. In most 
of all situations the combined interpolation yields a better approximation of 
A/f than one of the methods alone does. Furthermore, this method gives 
an unambiguous solution for the source interpolation problem. 

An even better alternative would be to compute the a, /3 and 7 coeffi- 
cients for a fixed combined linear and parabolic interpolation of the source 
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function, as in equation (4.5). The advantage is that the improved interpo- 
lation, as shown in figure 4.10, directly enters the analytic integral, instead 
of taking the mean of two inaccurate integrals. This has not yet been ac- 
complished in the scope of this work, but is planned for future work. 

4.6 Solving the rate matrix equation 

4.6.1 Computational weight 

As described in section 2.8.2 the rate matrix equation is a linear system 
of rank N, where N is the number of levels of all ionization stages of an 
atomic species. The number of levels per atom varies considerably between 
the species. In fact, the number is only limited by the atomic data that is 
used as input. The procedure for solving a linear system is straight forward, 
but for large systems this becomes computational expensive. 

The dynamic range in the solution of the system is very large ni € 
(0,Ne)- In praxis the dynamic range in the rate matrix is large too, and 
consequently the occupation numbers are sensitively coupled. For example, 
when determining rii, the transitions into level i from a weakly coupled, 
highly populated level j might be as important as the transitions from a 
strongly coupled, lowly populated level k. 

Examination of the solution of the rate matrix equations shows that the 
precision of FORTRAN 64-bit double precision variables is not high enough to 
obtain accurate results. For example, the insufficient precision causes the 
solution of the matrix equation (2.72) to contain negative n,. Therefore, 
in PHOENIX the rate matrix equation is solved using quad-double precision 
(QD). This is supported by an extension to FORTRAN-90, [HLB01]. A QD 
number is an unevaluated sum of four double precision numbers, capable of 
representing at least 212 bits of significand, or 64 digits. 

The higher precision causes large computational overhead. The memory 
required to store the rate matrix is a factor of 4 larger, which also degrades 
the CPUs cache performance (the cache hits/misses balance). Also, signif- 
icantly more floating point operations are needed per arithmetic operation. 
Added up, solving the rate equations is one of the most time consuming 
parts in the atmosphere calculation when large, realistic model atoms (with 
many energy levels) are used. 

4.6.2 Explicit rate matrix equation 

In order to discuss methods for solving the rate matrix equation (2.72) 
(RME) an explicit example of a rate matrix is shown in figure 4.11. The 
matrix corresponds to a simplified atom with four ionization stages. 

The matrix elements P are given by equation (2.70). These matrix 
elements connect states of different ionization stages. A colorized schematic 
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Figure 4.11: This is the explicit rate matrix for a simplified model atom, consisting of 
four ionization stages j. The neutral and doubly ionized stages (j = 1 and j = 3), apart 
from the ground states, have one excited state. The singly ionized stage has three excited 
states. The last stage is the bare nucleus. Within each ionization stage j the levels are 
sequentially numbered after their excitation energy i, starting from 1 for the ground level. 
The rates P are defined by equation (2.72). 

The black lines show the origin of the rates from the different ionization stages. In the 
squares the lines form around the diagonal, the rates originate from transitions between 
states of one specific ionization stage only. Rates between ionization stages only occur for 
the key levels that are ground state and continuum simultaneously. The key levels in this 
example are I2 and I3 (omitting the bare nucleus that has no further excitation levels). 
Note that the sum over each columns matrix elements is zero by construction. 



version of the explicit RME is shown in figure 4.12. This is the raw system 
that is not yet closed by the constraint equation (see section 2.8.2). Each 
ionization stage has its own color, and the matrix elements feature the colors 
of the ionization stages of the levels they are coupling. 



4.6.3 The old system 

In order to make solving the RMEs (one for each atomic species) computa- 
tionally treatable in PHOENIX a simplifying approximation is made. The rate 
matrix is almost block diagonal and can be made block diagonal with the 
following approximation. The ground states, which are at the same time the 
continuum state of the next lower ionization stage, are split up and treated 
in two parts: in one group as the ground state, and in another group as the 
continuum state. The block diagonal form of the rate matrix is 



/Pi 







p 2 







(4.8) 
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Figure 4.12: The explicit rate matrix of figure 4.11 shows matrix elements that couple 
levels of different ionization stages. In this schematic version of the RME each ionization 
stage j is coded with its own color: red, green, blue and red for j — 1, 2,3,4 respectively. 
The matrix elements show the colors of the ionization stages of the levels that they couple. 
Also, each m has the color of the ionization stage it belongs to. 

All diagonal elements have colors of multiple stages, yellow, turquoise, purple and gray, 
from an overlay of red-green, green-blue, blue-red and red-green-blue respectively. Also, 
the rows and columns corresponding to the key levels show mixed colors. The most mixing 
elements are the diagonal elements that correspond to the key levels. All other non-zero 
elements show only a single color. 

Note that this is the raw system that still contains one redundant equation (see section 
2.8.2). 



with j the last ionization stage of the atom. The off diagonal blocks are zero 
matrices. 

The RMEs fall apart into independent blocks, one for each ion of the 
atomic species under consideration. This system is shown schematically in 
figure 4.13 using the same color coding as in figure 4.12. Each block is a 
linear system for all levels of one ionization stage, plus the continuum of 
that ionization stage. Each system is closed by a constraint equation of the 
form 

c j c j 

i=l i=l 

with lj and Cj being the ground and continuum level of ionization stage j, 
and jVj the total number density of ion j plus continuum. In this notation 
the ground level of the neutral ionization stage li is the lowest energy level, 
denoted with 1 in equation (2.71). The system with constraint equations is 
shown schematically in figure 4.14. The constraint equations (4.10) provide 
the coupling between the separated ionization stages. 

The jVj do not directly depend on the rates. This allows to solve the 
linear system for each ion separately. The rank of each of the linear systems 
is much smaller than the rank of the full RME (2.72). Small systems are 
much faster to solve, which saves computation time. Furthermore, only 
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Figure 4.13: The rate matrix is almost block diagonal, see figure 4.11, and with a 
simplifying assumption it can be made block diagonal. The system falls apart into a 
small systems, one for each ionization stage. Those are easier to solve, deducing memory 
consumption and computation time. 

A great disadvantage of this method is that the simplification is not generally valid. The 
larger the deviations from LTE are, the more inaccurate the solution becomes. 
This equation is not yet of highest rank, it is closed by replacing redundant equations with 
constraint equations, see figure 4.14. 
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Figure 4.14: The block-diagonal RME of figure 4.13 is closed by constraint equations 
(4.10), one for each block. 
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relatively small amounts of memory are needed at a time, since the number 
of matrix positions scales with the square of the systems rank. 

Analysis showed that this method is not accurate for the physical con- 
ditions of typical atmospheres modeled in this work. The most obvious 
problem is that if departures from LTE are large, then also the ionization 
balance of the species under consideration strongly departs from the LTE 
balance. However, the ionization balance is not directly solved for. Itera- 
tively, the ionization balance can change, via the constraint equations. This 
is achieved by shifting the ion population from the ground state to the con- 
tinuum, or back, which is very ineffective if the initial ionization balance 
is far off from the true solution. Thus the convergence properties of this 
scheme are poor. The extra iterations required to roughly converge the so- 
lution of the RMEs are expensive and cost a multiple of the computation 
time gained by splitting the RMEs. 

The ionization stages are connected by the key levels that are the contin- 
uum for one stage and the ground state for the next higher stage. Splitting 
the rates of the key levels is physically only valid if the net rates to and 
from the key levels are exactly zero for each part of the rates separately (the 
parts related to the higher or the lower ionization stage). This requirement 
is satisfied by the solution of the RME, but that solution is just to be found. 

In the method artificial levels are added to the atom, the counterparts 
of the actual key levels. These additional levels require additional equations 
in order to close the system. Those are the equations (4.10). Mapping the 
results of this extended system back to the actual system is not possible, 
since the actual key level number density is not defined. Defining it as 
the sum or the mean of the twin levels generally violates number density 
conservation. Defining the key levels as just one of the pair, neglecting the 
other one is problematic, since there is no physical preference for one or the 
other. 

Only in the special case where the exact solution of the RME is used 
in the definition of the constraint equations (4.10), the number densities of 
both twins are exactly equal, and then the twins can be interpreted as one 
identical level. So when starting from the right solution, this method yields 
the right result. Especially, in the regions where LTE is a good approxima- 
tion the method gives reasonable results when starting with LTE as initial 
guess. On the contrary, where departures from LTE become larger, the ini- 
tial condition for the RME solver becomes worse, and thus the solution is 
not exact, and the rate of convergence degrades. 

This is shown by the example of a (typical) model of stage 3 (see table 
4.1) typical for this work in figures 4.15 and 4.16. In stage 3 Fe is newly 
added to the model and the population numbers rij of the other atomic 
species are kept fixed. The initial guess for the Fe population numbers is the 
Saha-Boltzmann distribution, which is unrealistic for the outer regions. The 
ionization balance from the first solution of the RME is plotted in figure 4.15, 
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Figure 4.15: This plot shows the Fe ionization balance of a model of stage 3 (see table 
4.1), in which Fe was just added. The m of all other atomic species are kept fixed to 
stabilize the radiation field while the rn of Fe adapt to it via the RME (solid curves), from 
their initial LTE values (dashed curves). The method to solve the RME that is used in 
this model is not exact. 

Except for Fe xviii in the outer regions (small layer numbers) all deviations from LTE are 
relatively small. 

In subsequent iterations the ionization balance is refined, see figure 4.16. 

in comparison with the initial LTE ionization balance. Except for Fexvin 
in the outer regions the departures from the LTE ionization balance are 
small. The change in the ionization balance in the following iterations, in 
which the rii of all other elements and the atmospheric structure are kept 
fixed, is shown in figure 4.16. The resulting ionization balance in the last 
of these 6 iterations is plotted in figure 4.17. Comparison of this iterated 
result with the exact solution from the next section, figure 4.19, shows that 
the overionization of the higher stages in the outer regions is slow, but the 
underionization of the lower stages is even slower. 

4.6.4 The full system 

The most straight forward, generally valid way of solving the linear system 
of rate equations is described in section 2.8.2. The schematic form of this 
system is shown in figure 4.18. One equation of the raw (unclosed) system is 
redundant and is replaced by the constraint equation (2.71). The choice of 
the equation to replace has a strong influence on the numerical precision of 
the solution. If, for example, the constraint equations replaces the equation 



4.6. SOLVING THE RATE MATRIX EQUATION 



61 
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60 
Layer 

Figure 4.16: This plot shows the weighted averaged changes to the Fe ionization balance 
AA^Fo, defined by equation (4.1), for 6 subsequent iterations. In each iteration the n» of 
all elements other than Fe are kept fixed, as well as the atmospheric structure. Clearly, 
in the outer regions, where the departure from the initial LTE ionization balance is large, 
the convergence rate is poor. The result from the subsequent corrections to the ri; in the 
inner regions is visible in figure 4.17, where the ionization balance after these 6 iterations 
is shown. 



62 CHAPTER 4. METHODS FOR GOOD NLTE MODELS 




120 100 80 60 40 20 

Layer 



Figure 4.17: The ionization balance of Fe after 6 iterations is quite different from the 
balance of the first solution of the RME. The LTE ionization balance (initial assumption) 
is shown by the dashed curves. Comparison with the exact solution of the next section, 
figure 4.19, shows that there are still significant differences. 

Around layer 35 the lower ionization stages Fexill-svi are by far not underabundant 
enough. In the inner regions, where the highest ionization stages Fe xxvi and Fe xxvn are 
yet out of LTE, the whole ionization balance is driven out off the correct LTE balance. 
Because of the poor convergence properties and the extra iterations needed with this 
method this method is replaced by a new one, described in section 4.6.5. 
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Figure 4.18: The system of rate equations is closed by the constraint equation of number 
density conservation. This equation replaces one of the redundant rate equations. 
The choice of equation to replace has an important influence on the numerical accuracy 
of the solution. A good choice is a key level k with a large expected in the solution of 
the system. 
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Figure 4.19: This Fe ionization balance (solid curves) is computed with the full RME 
(4.12), for the same model as in figure 4.15. For comparison also the LTE balance is shown 
(dashed curves). Solving the full RME is computationally very expensive, but yields the 
exact solution. This solution is used to check the results from other, more resource saving 
methods. 



of matrix row k and the number density is very small <C Ne, then 
analytically 

n k = N E -J2 n i ( 4 - 13 ) 

is the difference of two almost equal numbers. This choice is very prone to 
errors from the limited significance of numerical variables. 

Experiment shows, that the best choice is the key level with the largest 
number density. This is not known prior to solving the system but a guess 
can be based upon a Saha-Boltzmann distribution, or better, if available, 
on the solution of the system from previous iterations (the current n^). 

This method yields the exact solution of the RME. The resulting ion- 
ization balance from the first solution of the RME for the same model as in 
figure 4.15 is shown in figure 4.19. The departures from the original LTE 
ionization balance are much larger than with the old method of the previ- 
ous section. In subsequent iterations, in which the nj of other elements and 
the atmospheric structure are kept fixed, the Fe ionization balance hardly 
changes. 
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Figure 4.20: After the first iteration the weighted averaged changes to the Fe ionization 
balance AA^Fe are very small in all regions of the atmosphere. It shows that despite strong 
deviations of the rii for Fe, the radiation field was good enough to yield accurate rates for 
the largest m (that weigh most in the average). This results from keeping all other m 
fixed when adding Fe. 

4.6.5 The new system 

For the reasons discussed in section 4.6.1 the straight forward method of the 
previous section is computationally very expensive. A different approach 
is required in order to reduce memory consumption and calculation time 
to a level appropriate for presently available computer power. The follow- 
ing method provides both an exact solution of the RME and a significant 
reduction in memory consumption and CPU time. 

The raw RME (2.72) without the closing constraint equation (2.71) does 
not have a unique solution. The solution of the closed system n satisfies the 
open system, but also the scaled solution an does, with a some constant. 
So if any such vector aft can be found, then it can be scaled to the unique 
solution n using the constraint equation. Thus a different constraint can be 
used to solve the matrix. 

One possibility is the simple constraint of fixing one of the number densi- 
ties to some constant = 1. If the level k that is being fixed is chosen to be 
one of the key levels, then this completely removes the dependency between 
the ionization stages above and below that key level. In other words, this 
choice splits the linear system in two independent parts. A block diagonal 
matrix is obtained, with two blocks. This method can be called the 2-block 
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Figure 4.21: The unclosed system, figure 4.12, does not have a unique solution. However, 
since its right hand side is the zero vector b — 0, its linear properties allow to quickly find 
one arbitrary solution. That arbitrary solution of the unclosed system can then be scaled 
to satisfy the closing constraint equation. With this method the unique solution is found 
in a much faster way than when the closed system is solved. 

The modification to the RME that quickly yields one arbitrary solution is shown here. 
One redundant equation is replaced by a fixation of one of the key variables. Thus the 
system of rank 9 becomes block diagonal, and splits into a system of rank 6 and a system 
of rank 2, which are solved faster than the system of rank 9. 

Large atoms with many ionization stages possess many key levels, which brings the ranks 
of the two blocks closer together. So the splitting is more efficient for large atoms than 
for small atoms, like in this example. However, those are already easy to solve. 



diagonal grouping method. This is shown schematically in figure 4.21. 

The maximum gain from this method is achieved if the system can be 
split in two parts with equal rank. In that best case scenario the memory 
savings are a factor of 4. In praxis the most centered key level is often not 
exactly in the middle, especially for small atoms, with only few ionization 
stages. But small atoms are easy to solve. More importantly, for large atoms 
with many ionization stages, the gain is often not far from the theoretical 
maximum. 

The savings in computation time needed to solve the system are even 
better. Comparative tests performed for Fe, by far the largest atom used 
in this work, showed that the CPU time with the new method is a factor 
of 7-15 smaller than with the full system method of the previous section. 
The exact gain depends on the ionization balance in the converged solution 
of the RME. The highest gain is achieved in the case where the ionization 
balance is centered right at the centermost key level m\ that was kept fixed, 
i.e. if in the solution n mi > n\ for any other key level I ^ m\. 

This allows to fine tune the choice of key level to fix. The exact center 
lies (most often) between two key levels. From the current population 
the best populated key level I in the solution of the RME is guessed. If the 
second best centered key level wi2 is not much further off center than the best 
centered key level mi and m,2 is closer to I than mi, then mi is kept fixed 
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rather than m\. This yields a slightly worse gain in memory consumption, 
but a significant better gain in required CPU time. 

Because the CPU time depends on the ionization balance it depends on 
the region in the atmosphere. Therefore, if the RME solver is parallelized 
over layer, then a round robin scheme is better than a clustered (blockwise) 
spreading. Thus the layers with a disadvantageous ionization balance are 
distributed over the workers. 

4.6.6 Dealing with numerical inaccuracies 

As described in section 4.6.1 the dynamical range in the rate matrix and 
in the level populations is large. The RME are first solved using the LU 
factorization method. This is an analytical method. The accuracy is only 
limited by the numerical precision of the variables with a given condition 
number of the matrix 1 . Iterative linear system solver methods, like the 
Jacobi method or the Gauss-Seidel method start from an initial guess that 
is iteratively refined. These methods are very slow 2 for typical large RME if 
the initial guess is poor. Therefore, they are not suitable to solve the RME 
from scratch. But they are very useful to refine the result from the analytic 
solver. The numerical inaccuracies (like negative rij) are removed from the 
solution of the RME within a few steps. 

Thus the combination of the analytical solver and an iterative solver 
is very powerful. And another benefit can be drawn from the solver com- 
bination. Since numerical inaccuracies from the analytical solver can be 
corrected afterwards, this allows to reduce the precision of the variables 
used by the analytical solver. Instead of QD variables double-double preci- 
sion (DD) or even standard double precision (SD) can be used. Afterwards 
the precision is checked with the iterative solver, and if the inaccuracies are 
too large to be corrected within a few iterations, then the system is solved 
analytically from scratch with a higher precision. The computational over- 
head of higher precision variables is so large, that the wasted CPU time of 
the lower precision solver step is relatively small. Going from QD to DD 
saves a factor of 7 in CPU time, and going from DD to SD saves another 
factor of 20-30! The savings can be attributed to less flops per arithmetic 
operation, better cache performance, and better performing compiler opti- 
mization (especially when going from DD to SD). These savings combined 

1 The condition number is a property of the matrix itself, not of the algorithm or the 
floating point accuracy. It specifies how numerically well-conditioned the equation is. In 
the case of the linear system Ax — b, a high condition number (ill-suited for numerical 
computation) would mean that even a small error in b would yield a large error in the 
solution x. 

2 When a Matrix is block diagonal the solution of each block of the matrix equation is 
independent of the other blocks. In the RME considered here, the matrix is nearly block 
diagonal, so that the solutions of the blocks are not fully independent, but the coupling 
between the blocks is yet (very) weak. 
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Figure 4.22: An example of oscillations that occur in realistic models when an undamped 
UL temperature correction is used. Here the relative correction to the Temperature is 
plotted for 64 atmospheric layers (spanning the optical depth range 10~ 10 < r < 10 2 ) 
over 10 model iterations. 

with the gain from the 2-block diagonalization method of section 4.6.5 result 
in an extreme reduction of CPU time. 

With this result, the CPU time consumption for solving the RMEs, in- 
stead of being the dominating and limiting part of the atmosphere model 
calculation, becomes almost negligible against the other parts (e.g. opac- 
ity computation and radiation transport). This result is also important 
for other massive NLTE models, like supernova models, where the highly 
abundant iron-peak elements Fe, Co and Ni have very many levels. 

4.7 The Active Damping method 

One of the challenges of atmosphere modeling is to find a correction proce- 
dure that converges for a large variety of atmospheres. Normally, the price 
for such universality is that the rate of convergence is not optimal for dif- 
ferent types of atmospheres. But a good convergence rate is important to 
save calculation time. 

Analyzing the Unsold-Lucy (UL) temperature corrections for realistic 
models shows that undercorrections and overcorrections (oscillations) occur. 
An example of oscillation is given in figure 4.22. A number of damping and 
acceleration methods were proposed in the literature, e.g., [Luc64], [Dre03], 
[HG03] , but all of them either need to be individually tuned or are ineffective. 
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Figure 4.23: Three cases of damped harmonic oscillators: 1. Underdamped, oscillation 
occurs; 2. Overdamped, the equilibrium is approached slowly from one side; 3. Critical 
damping, the equilibrium is approached in the fastest possible way. 

In this section a new "Active Damping" (AD) method is described, in which 
the optimal damping factor for each layer and each iteration is determined 
on-the-ny from the values of previous corrections. 

In classical mechanics, a system with a restoring force can be over- 
damped, underdamped or critically damped, see figure 4.23. The UL cor- 
rection acts like a restoring force to the temperature driving it towards 
equilibrium and in analogy there exists an optimal damping factor to this 
correction. 

4.7.1 Active damping for linear corrections 

As a first approximation the UL correction AB, in the following written as 
A, is assumed to be positively proportional to the deviation d from the final 
temperature (for which energy conservation is satisfied). 



with / the amplification factor and the index i labels the iteration number. 
In the case / = 1 the undamped UL correction would converge within a 
single iteration. In praxis, this rarely happens to single layers, let alone to 
the complete atmospheric structure. 

In order to develop the formalism, assume that a damped correction A^ 
is applied 



A = / • <k 



(4.15) 



Pi >0 



(4.16) 
(4.17) 
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with ^ being the damping factor for iteration i. The restriction to positive 
/3 reflects the assumption that A is positively proportional to the deviation 

d. After the first correction d becomes 

d 2 = d 1 + A[ = (l-f3 1 f)d 1 (4.18) 
and using equation (4.15) the second correction becomes 

A 2 = -/(l - frf)di (4.19) 
For / one finds the expression 

Ai - Ao Ai - A 2 , 

The AD-factor a is defined to compensate for deviations from / = 1, so that 
the correction damped with a cancels d completely 

ai Ai = -di (4.21) 

In order to calculate the AD-factor, values of two subsequent UL corrections 
are used. Using equation (4.15) and (4.20) and generalizing for all iterations 
one finds 

where is the generalization of equation (4.20) 

fi = ^ (4-23) 

Pi-l^i-l 

The damping factor can now be set equal to this AD-factor aj. In the 
first iteration, as Aj_i is not yet known, some preset value must be used, 

e. g. (3 = 1 for an undamped correction. Note that in equation (4.15) 
/ is assumed to be constant, whereas in equation (4.22) /j changes from 
iteration to iteration, /j is the constant value that fits to iterations i — 1 and 
i. With this value of /, the UL corrections for these two iterations would 
be proportional to di, as in equation (4.15). This proportionality is never 
exact, so that after the second iteration the UL correction is found to be 
not yet completely converged, thus A3 / 0. 

4.7.2 Active Damping for non-linear corrections 

The assumption that the UL correction depends linear, equation (4.15), can 
be rewritten as 

Ai = f-s d \di\ (4.24) 
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with 



s d = sign(dj 



(4.25) 



Assume that the UL correction satisfies a power law instead of a linear 
relation 



Aj = / ■ s d \di 



(4.26) 



where < 7 < 00. Note that dk and its sign are not known a priori. In 
analogy to the linear case, equation (4.16), in every iteration the correction 
Aj is damped with some factor 



A ■ = A ■ /j ■ s d \di 



(4.27) 



fi and 7j are the values of / and 7 for the current iteration i. Multiple itera- 
tions will be needed since, like in the linear case, the non-linear assumption 
does not hold exactly. Some algebra yields the expression for fi 



Si-ilAi-ilVT-sJAJVT 



with 



Si-! = sign(Ai_i) 
Si = sign(Ai) 

' s i -!\A i _ 1 \ 1 /'r - Si\Ai\^ 



Sf = sign 



A-lA;_i 



(4.28) 



(4.29) 
(4.30) 

(4.31) 



The parameter 7, is determined by the requirement that (in the assumption 
of equation (4.26)) the amplification factor f(j) is a constant for some value 
of 7 



/i-i(7) = Ml) 



(4.32) 



This equation can be solved numerically using the UL corrections of three 
subsequent iterations Aj_2, \-i arid Aj, which yields both fi and 7, in the 
following called 7^. 

Finally, the expression for the AD-factor a in the non-linear case is 



OLi. = 



-dj 
A, ; 



(4.33) 



Note that for 7, = 1 both fi and ccj reduce to the linear versions of equations 
(4.23) and (4.22). 
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4.7.3 Implementation 

When the AD method is applied to stellar atmosphere models a few difficul- 
ties arise. These are caused by the fact that the assumptions of equations 
(4.15) and (4.26) do not hold exactly. In fact, there exists no relation that 
holds exactly because of the assumptions of the Unsold-Lucy scheme. 



Problems 



One problem is related to the singularity in equations (4.22) and (4.33) 
for Aj = Aj_i. The AD-factor behaves asymptotically in the neigh- 
borhood of the singularity. Asymptotic jumps often occur when the 
corrections A become small, as in this case numerical inaccuracies in 
A become large. These jumps in the correction factor have no physical 
justification. 

Negative AD-factors occur when the UL correction gets larger than in 
the preceding iteration and is in the same direction (either positive or 
negative) 



sign(A^ 
sign (A 



i-i 



| Ai| > |Ai_i| (4.34) 



Assuming that the UL corrections are always restoring, i.e. directed 
towards the equilibrium, these negative corrections are invalid. 

3. Large jumps in the AD-factors of adjacent layers cause unphysical 
jumps in the resulting temperature structure. 

4. Large AD-factors (large extrapolations) yield inaccurate results and for 
very small factors the calculation of the subsequent AD-factor becomes 
inaccurate. 

5. The solutions 7 of equation (4.32) can be inappropriate, i.e. 7 > 1 or 
7<C 1. The UL correction is not extremely non-linear. 



Basic solution scheme 

The following scheme proved to handle the problems described above: 

1. Compute for all layers /. 

2. Reject extreme values < 0.5 or 7 (I) > 1.3. 

3. Treat layers without valid 7 as linear = 1. 

4. Smooth 7 according to -y(l) = - 1) + -y(l) + j(l + 1)) /3 

5. Compute /(/) using the smoothed 
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6. Mark 5 layers 3 (I — 2, ■ ■ ■ ,1 + 2) around each layer where f(l) < 0. 

7. Linearly interpolate / over the marked ranges. 

8. Calculate the AD-factors a(l). 

9. Limit the new damping factors (3(1) according to 0.25 < (3(1) < 2.5 

10. Finally the new damping factors (3i(l) should be smoothed over 5 layers 
according to 

m _ 4^+2^-l)+2 W + l) + W -2) + W + 2) (4 35) 

All limiting values in this scheme were determined experimentally. 
Limitations 

In some cases the basic scheme fails 

1. step 7 is not possible for ranges that extend to the outermost /innermost 
layer. 

2. if the range to be interpolated over in step 7 is very large, i.e. larger 
than #layers/4 

3. if layers with large jumps of the resulting damping factor (3 exist, i.e. 
when two or more of the following four ranges are exceeded 4 : 

2/3 < \p(l)/P(l±l)\ < 3/2 (4.36) 
l/2<|/3(/)//3(/±2)|<2/l (4.37) 

Instead of AD the following more primitive method can be used to calculate 
the damping factor. 

4.7.4 The alternative: Passive Damping 

In the occasions that active damping fails, one has to revert to another 
method to determine an appropriate damping factor. The primitive method 
described in the following is called "passive damping" (PD), indicating that 
the damping factor is not intended to completely cancel the deviations d in 
one step in contrast to active damping. The scheme for passive damping is 
as follows 



3 When less than 100 layers are used, it may be better to mark 3 layers around layer I. 

4 This is appropriate for models with up to 100 layers. When more layers are used, 
values closer to unity (e.g. 4/3 and 3/2 respectively and their inverses) may give better 
results. 
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1. For each of those layers I the damping factor of the previous iteration 

is taken, and discretized to the values [0.3, 0.42, 0.6, 0.8, 1.0, 
1.25, 1.6, 2.0]. 

- If strong oscillation occurs (—1 < Aj/Aj_i <C 0), the value is stepped 
down by one. 

- If strong overdamping occurs (0 <C Aj/Aj_i < 1), the value is stepped 
up by one. 

- If oscillation or overdamping is very small (|Aj/Aj_i| <C 1), the value 
is retained. 

- If the correction increases (|Aj/Aj_i| > 1), interpolate the damping 
factor (on the discrete grid) over this layer using the two closest layers 
for which this is not the case. For the outermost /innermost layer then 
the center value of 0.8 can be used. 

2. Now the differences in the damping factors between the layers must 
be reduced to be one step on the discrete grid at most. This is done 
starting from the two enclosing layers for which the AD-factors are 
still used, going towards the center of the range. Only in case there 
are less layers than steps needed to go from one border to the other, 
differences larger than one step are to be left over. 

4.7.5 Convergence properties 

In the calculation of model atmospheres the UL temperature correction is 
found to converge much faster with AD than with fixed damping constants. 
For all types of models tested thus far (white dwarf, brown dwarf, red giant, 
nova and supernova models) a speed up of a factor of two or more is achieved. 

Very quick models, i.e. models that converge within 1-3 iterations with- 
out the use of AD, do not benefit from the new method, because it needs 
1-2 iterations to initialize (for the linear or the non-linear case respectively). 

For the comparison of convergence rates, data points in two dimensions 
must be analyzed, one dimension for spatial layers and one for temperature 
correction iterations. A convergence of all model layers is required. However, 
an error in a single layer is less significant than an error extending over 
multiple layers. 

In figure 4.24 the convergence properties are shown for a hydrostatic 
LTE model. This type of model is suitable for the purpose of demonstra- 
tion. Firstly, in a hydrostatic atmosphere, being relatively compact, the 
diversity in physical condition in the atmosphere is smaller than for the ra- 
dially extended type atmospheres introduced in chapter 3. Furthermore, the 
NLTE generalization of the Unsold-Lucy temperature correction method, 
see section 4.8.1, requires special care to deal with the outer regions in the 
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Figure 4.24: These plots show typical model convergence properties of the UL correction 
with different damping methods: Active Damping (top), undamped (/3(l) = 1) (middle) 
and constant damping /3(l) = 0.5 (bottom). On the left the undamped relative correction 
A oc AT/T computed with the Unsold-Lucy correction procedure is shown for 64 model 
layers and 10 iterations. On the right the same data is plotted logarithmically log |A| on 
a scale from 10" 5 to 10 _1 . 

The rate of convergence with AD is larger than with conventional constant damping: after 
3 iterations an accuracy of 0.6% is reached, for which 9 iterations are needed in case of 
j3 — 0.5 and even more in the undamped case. The large rate of convergence is maintained 
even when the corrections become very small. After 10 iterations the largest A with AD 
is a factor of 100 smaller than the largest A without AD. 
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atmosphere, where the coupling between the radiation field and the thermal 
pool is very weak. These procedures have been implemented in the active 
damping part of the code, and not in the classical globally constant damping 
methods. Therefore, the models typically do not converge at all with the 
constant damping methods. But the effect achieved with active damping, as 
shown by example in figure 4.24, is found in every type of stellar atmosphere 
model. 

The sample model in figure 4.24 has been calculated for three cases: using 
a damping factor of (3(1) = 0.5, undamped (3(1) = 1. and actively damped 
(3(1) = a(l). (3 = 0.5 has proven to be a good value, being small enough to 
cancel most oscillations. For simplicity, let us call a model converged when 
the corrections get below 1% in every layer (for production models more 
detailed convergence criteria are used). Then 8 iterations are needed for 
(3 = 0.5, 18 iterations for (3 = 1 and 3 iterations with AD. So in this model 
the convergence speed up is about a factor of 2.5. 

Note that (3 = 0.5 is not fine-tuned, a factor of (3 = 0.7 yields convergence 
within 5 iterations for this simple model. 

4.7.6 Conclusions 

The Unsold-Lucy temperature correction can be improved significantly by 
appropriate damping. A constant damping factor can be fine-tuned in order 
to optimize the rate of convergence. However, as the situation changes from 
iteration to iteration, a constant damping factor never is optimal. Further- 
more, different regions in atmospheres behave differently, which is again not 
supported by a constant damping factor. 

For the majority of atmospheres a speed-up of about a factor of two or 
more is achieved with the Active Damping method, without the need for 
any fine tuning. However, for simple models that converge within a few 
iterations using a constant damping factor, the gain is lower than the factor 
of two, because the method needs 1-2 iterations to initialize. 

Generally, the Active Damping method is not limited to the Unsold-Lucy 
temperature correction, but can be applied to any physical system in which 
a single independent quantity must be brought into some equilibrium. 

4.8 NLTE generalization of the Unsold-Lucy tem- 
perature correction method 

4.8.1 The formalism 

The derivation of the Unsold-Lucy method (section 2.11.1) requires the rep- 
resentation of the source function S as a 'weighted' sum of B and J (equation 
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(2.105)) 

S x = (k x B x + a x J x )/(k x + a x ) = eB x + (1 - e) J x (4.38) 

e = KX (4.39) 
+ e\ 

This expression is valid for LTE, but not for NLTE. For example, in NLTE 
the source function S can exceed the sum of B and J (and thus also the 
weighted sum of equation (4.38)) or fall below both B and J. Clearly, this is 
not compatible with the assumed expression of the source function, equation 
(2.105). 

In general the (macroscopic) source function is built from the (micro- 
scopic) rates. The rates are functions of the radiation field, the level pop- 
ulations, and the electron density, and can not be attributed to thermal or 
non-thermal (see section 2.3.3) processes. Hence S cannot be represented in 
the classical form, equation (2.105), if the rates are used to calculate S. 

However, if, in the classical view, at a specific location in the atmo- 
sphere thermal processes play a dominant role over non-thermal processes, 
then the atmosphere is forced into a state close to LTE. On the other hand, 
if non-thermal processes are dominant, then the departures from LTE will 
become large. Now, reversing these implications, one can estimate the bal- 
ance between thermal and non-thermal emission by the order of departure 
from LTE. 

Writing the emissivity as the sum of a thermal and a non-thermal fraction 

Vx = V? + vT h (4-40) 

allows to use equation (2.106) as it stands, but now with the averaged opac- 
ities kb and kj written as 

"J = jJ XxJx~vT h dX (4.42) 

and xh is defined by equation (2.110). 

Let each bound-bound transition have a statistical fraction 9 of emission 
that is thermal 

bb _ bb,th . bb.Nth 

Vx = V\' +V X 

^(«) (4.43, 
= Yl + (! - 
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where i and j are the levels of all transitions that contribute to the emissivity 
at wavelength A. Then the thermal and non-thermal terms of equation (4.40) 
can be written as 

Vf = 4 + ^ + ^' th (4.44) 
= axJx + ^Nth (445) 

The emissivity from bound-free and free-free transitions can be assumed 
to be fully thermal. In both cases the interaction takes place with a free 
electron that has a Maxwellian velocity distribution [Mih78]. Note that 
the NLTE expression for free-bound emissivity (2.68) depends on the LTE 
occupation number of the bound state to which the recombination takes 
place. 

For bound-bound transitions the thermal fraction Oij can be estimated 

as 

9 l3 = Twn(bi/bj,bj/bi) = mm(b*/b*,b*/b*) (4.46) 

The minimum function reflects the fact that departures from LTE can man- 
ifest as underpopulation or overpopulation. If the upper level and the lower 
level have the same departure b from LTE, then for the transition a pseudo 
LTE situation exists. The thermal fraction must then be close to unity. 



4.8.2 Inferior alternatives 

The new NLTE generalization of section 4.8.1 extends the applicability of 
the Unsold-Lucy temperature method to situations where the emissivity is 
not fully thermal. 

Originally, in PHOENIX there were two distinct ways to treat the NLTE 
opacities in the temperature correction procedure. In the first variant the 
NLTE opacities are left out from the integrals of the expressions for kb, 
kj and xh (equations (2.108) to (2.110)). Since the non-thermal opacities 
originate from NLTE opacities only, this circumvents the problem. 

This is a useful approximation if, as usual, the bulk of the opacities is 
LTE opacities. But, as described in section 4.2, the models performed in 
this work are pure NLTE models, meaning that there are no LTE opacities 
at all. However, the non-thermal opacities originate from the bound-bound 
transitions only. Therefore, in analogy to the original idea, the NLTE bound- 
bound opacities are ignored, and the free-free and bound-free opacities are 
used as approximation in the opacity means. From equations (2.110), (4.41), 
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(4.42), (2.25), (4.44) and (4.45), omitting the bound-bound terms, it follows 
KB = bS ^ + ^ dA= l / xfB x + V ^dX (4.47) 
Kj = lj (^ + X^ f )«/AdA (4.48) 
) H x d\ (4.49) 



\// = -J i h x -XA b 



In the second original variant the NLTE emissivities are treated like 
LTE opacities so that the opacity means Kb and kj (equations (2.108) and 
(2.109)) become 

KB= B J(Xx-^)B x d\ (4.50) 
«J = -j j (XA-^A)JAdA (4.51) 
where % is given by equation (2.25). xh is defined by equation (2.110). 



4.8.3 Comparison 

In order to test the assumptions made in the NLTE generalization of section 
4.8.1 the method must be tested and compared to the alternatives of section 
4.8.2. In all variants the temperature corrections are determined by the 
same equations (2.120) and (2.121), but different definitions for the opacity 
means are used. 

The temperature structures obtained from the three variants are plotted 
in figure 4.25. In the optically thin outer regions the temperature struc- 
tures differ significantly. But in the deeper ranges they are almost identical. 
There the departures from LTE are small, and thus the non-thermal part 
of the emission goes to zero. Note that in LTE, where all hi = 1, the NLTE 
generalization of section 4.8.1 reduces to the second original variant. Fur- 
thermore, J\ becomes similar to B\, so that the definitions of kb, kj and 
Xh become identical. 

One important test is if the temperature correction method converges 
to radiative equilibrium, as described in section 2.10. For the model shown 
in figure 4.25 the two measures for relative deviations from RE (defined by 
equations (2.100) and (2.101)) are plotted in figure 4.26. After a number 
of temperature correction iterations the errors converge. The errors with 
the two alternative methods are significantly larger than those with the new 
NLTE generalization method of section 4.8.1. 



4.8.4 Limitations of the NLTE generalization 

The NLTE generalization of the Unsold-Lucy temperature correction method 
treats the opacities in a physically better way than the old alternatives. It 
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Figure 4.25: The temperature structures, plotted against the density, for a typical NLTE 
model resulting from three different variants of temperature correction procedure: 1. the 
new NLTE generalization of section 4.8.1 (black), 2. an approximation neglecting all 
bound-bound opacities (red), and 3. an approximation where the emissivities are treated 
like in LTE (green). 

In the inner regions, where departures from LTE are small, the three methods give very 
similar results. In the outer regions departures are large but with the small densities 
the opacities are small. The most important differences are in the middle regions, with 
densities between 10 -4 and 10 -10 g/cm -2 . 



was shown that also the convergence performance is superior. These two 
properties make the NLTE generalization the preferred temperature correc- 
tion method, and throughout this work this method is used for all model 
computations shown in this work, unless specified otherwise. 

There are, however, still a few limitations that are inherent to the UL 
method. In the UL method it is assumed that the proper temperature of 
the gas can be derived from the radiation transport equation (see equa- 
tions (2.102) and the following). This can only work if the opacities do 
significantly dependent on the temperature. As shown in section 4.4 there 
are atmospherical conditions in which such temperature dependence is not 
given. This is especially the case for regions where the departures from 
LTE are large. The population numbers are then so far off from their TE 
values, and so are the opacities, that it becomes very difficult to derive the 
temperature from it. 

In the NLTE generalization scheme this problem is apparent in the def- 
inition of kj, equation (4.42), where the integrand can become very small 
or even negative if the non-thermal part of the emissivity becomes large. A 
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1 2 3 4 5 6 120 100 80 60 40 20 

Iterations Layer 

Figure 4.26: These plots show the radiative equilibrium errors for the same model 
shown in figure 4.25. The errors are defined in equations (2.100) and (2.101) as £RE,int, 
corresponds to the integral form (H — Ho), and £RE >( uf, corresponding to the differential 
form (S — J) of the RE condition. The left plot shows the root-mean-square (RMS) 
of the errors over all model layers for subsequent temperature correction iterations on 
a logarithmic scale. The errors in the last iteration (number 6) are plotted against the 
model layer number on the right. 

After a number of iterations the errors converge. The errors with the two alternative 
methods are significantly larger than those with the new NLTE generalization method of 
section 4.8.1. The layers 40-100 have the most important influence on the spectrum, being 
between optically thick and thin. Right in this range the NLTE generalization method 
performs better. 
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small difference of two large numbers is very prone to errors. But the UL 
method is based on the assumption that the ratio of the opacity integrals, 
kj/kb, are independent of the temperature, i.e. constant over temperature 
correction iterations. In such cases, where the error in kj becomes large the 
formalism cannot work accurately. 

However, in those cases where the UL method fails, because the opacities 
(not the integrals) are insensitive to the temperature, the temperature is 
unimportant for the state of the gas, and so for the radiation transport. 
Therefore, where the temperature cannot be derived accurately, there the 
temperature does not need to be determined accurately. 

A 'resort' temperature method 

In order to determine the temperature approximately in those regions where 
the UL method fails, another method must be resorted to. A very simple as- 
sumption can be made using the Stefan-Boltzmann law (equation (2.94)) for 
a luminosity that is constantly extrapolated for outer regions. The following 
temperature relation follows 



with r re f being the outer reference layer for which the UL scheme is still 
working reliably. In temperature structures shown in this work this is ap- 
parent by a sudden very smooth, almost linear decline of the temperature 
against the density in double logarithmic plots. For example, in figure 4.25 
this occurs from layer 6 outwards, where layer 6 was the reference layer (the 
outermost layer that was yet determined from the NLTE generalized UL 
method). 

Validation condition 

One thing left to define is a validation condition for the UL method. The 
following condition has proven to work well. The UL method can be applied 
if the opacity integrals are approximately constant. This is tested by the 
condition that the opacity integrals kj and xh need to have been reasonably 
constant for the actual and the previous temperature correction iteration, 
where for reasonably constant these must satisfy both 
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The weaker condition for kj reflects the fact that the uncertainties in kj are 
larger due to the difference in its definition, equation (4.42), as discussed 
above. 

4.8.5 Line cooling/heating effects 

The monochromatic integrand of the radiative equilibrium condition, equa- 
tion (2.98) 

x\ = Vx ~ X\J\ = X\{S\ - JA) (4.56) 

here denoted with x, is called the net radiative cooling rate, see also [Rut95]. 
When S\ > J\, then locally an overdose of radiation is produced, repre- 
senting an energy loss for this wavelength. If x integrated over a line profile 
is positive then the line has a local cooling effect, called line cooling. The 
opposite is called line heating. Also, the continuum can have a local cooling 
or heating influence. 

In figures 4.27 to 4.29 this line cooling effect is demonstrated. If line 
cooling/heating gets dominant over the contribution of continuum cool- 
ing/heating then deducing the temperature from the interaction with the 
radiation field becomes difficult, since the line cooling/heating contributions 
strongly depend on the run of the mean intensity over the line profile. In 
the case of NLTE, the radiation field is iteratively brought into equilibrium 
with the radiative rates so that especially in the line centers the mean in- 
tensity is subject to changes. This further complicates the derivation of the 
temperature from line cooling/heating effects in regions that are optically 
thin apart from the lines that contribute to the cooling/heating. 
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Figure 4.27: These graphs show the effect of line cooling (see also the next figures). 
In the upper graph the gradual development over wavelength of relative radiative equi- 
librium error £RE,dif (equation (2.100)) is plotted. The nominator and denominator are 
plotted cumulatively (|cum) (red and green). The ratio of these sums is evaluated per 
wavelength (black). 

The lower graph shows something similar to the upper graph. Here not the RE is evalu- 
ated, but the UL temperature correction, specifically only the ABi term, equation (2.120), 
that is important in the outer regions where line cooling effects come into play. The con- 
stituents B ■ kb and J • Kj are plotted cumulatively (red and green). The ratio of these 
two minus I (equal to ABi/B) is evaluated per wavelength (black). 

The middle graph shows the integrands of ftp and kj (not cumulatively) on a logarithmic 
scale (red and green). Negative values of the latter are plotted positively in blue. 
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Figure 4.28: This figure shows a zoom in from figure 4.27 to a smaller wavelength range. 
The optical depth for this layer is r rc f = 2 • 10 -4 . At this depth, the atmosphere starts to 
become optically thin. 

Comparing the upper and lower graphs of figure 4.27, it is seen that the relative RE error is 
more sensitive to line cooling effects than the UL temperature correction. For this reason 
the UL temperature correction is more stable and features better convergence properties 
than a Newton- Raphson scheme based on RE. 

The lower graph shows that the strong lines in this wavelength range are heating whereas 
the continuum is cooling. Note that in the upper graph, the continuum is also cooling 
(the RE err increases) but here the lines are cooling instead of heating. 
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Figure 4.29: Here the same plots are shown as in figure 4.27, but now for a layer 
further out, with optical depth r rc f = 2 ■ 10 _J . At this shallow depth the lines start to 
dominate the heating/cooling influence of the radiation field to the matter. In this region 
also the UL temperature correction starts to become unreliable, because of inaccuracies 
in the mean intensity in the line centers of the lines that dominate the cooling effect. 
The problem arises at wavelengths where the monochromatic contribution to kj becomes 
strongly negative. This effect is inherent to a weak coupling of the radiation field to the 
matter in NLTE models. It is discussed in section 4.8.4. 
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Chapter 5 

First results 



The primary accomplishment of this work is the framework described in 
the previous sections. This framework allows to compute radially extended, 
massive (atoms levels) NLTE models. Once this framework is ready and 
working, it remains to show that it produces useful results. For this pur- 
pose a modest amount of models has been computed, as many as available 
computational power allowed for in the time frame of this work. These first 
results are shown in this chapter. 

An important test for atmosphere models is if systematic variations of 
single model parameters produce systematic results. This test was per- 
formed with a few small model series (model grids) . The size of these grids 
is yet small, limited by computation time, and further expansion is planned 
for near future. Apart from the purpose of testing the new framework, these 
grids give the opportunity to explore the parameter space of the models. 
Furthermore, careful analyzation of the impact of single model parameters 
on the results is an important preparation for a detailed spectral analysis of 
the available observational data, for which the models are fine-tuned to each 
observation individually. 

In section 5.1 the results are shown for a few small model grids, for 
both the hydrostatic and the hybrid type (nova) atmosphere models. Sub- 
sequently, in section 5.2 LTE and NLTE models are compared in order to 
confirm that a NLTE treatment is necessary for this type of models. Finally, 
in section 5.4 the models are compared to the currently available 'high-res' 
X-ray observations of novae in the super-soft X-ray state (V4743 Sgr, RS 
Oph, and V2491 Cyg). It is stressed at this point, that so far only a small 
amount of models has been computed. Therefore, the models are far from 
fine-tuned to fit the observations. Nevertheless, these first results fit the 
observations remarkably well, and a vast improvement is achieved with re- 
spect to previous work with PHOENIX [Pet05], where models were fine-tuned. 
These results can be seen as a very promising starting point for a thorough 
analysis of the observations through a detailed modeling process. 
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5.1 Model grids 

So far, parameter grids have been computed around the most basic model 
atmosphere parameters. 

Grids are computed for pure hydrostatic models and for the hybrid mod- 
els as they are described in chapter 3. The former type of models has two 
advantages. Firstly, there are only two basic model parameters (T e g and 
log(g)) in contrast to five for the hybrid models (T c g, log(g), M, Voo and f5). 
Therefore, the dimension of the grid is much smaller and thus the number 
of models to be computed. Secondly, the computational demands are signif- 
icantly lower than for the nova-type hybrid models. The two main reasons 
for this are: 

1. Hydrostatic atmospheres are relatively (radially) compact, so there is 
less diversity in physical conditions in the atmosphere. Consequently, 
a smaller number of layers is needed to sample the atmosphere. 

2. Since the velocity field is zero there is no coupling between the wave- 
length points so that numerically each point can be treated indepen- 
dently. This is an ideal situation for parallel computing. It features 
almost perfect scalability (to multiple parallel tasks). 

For the hydrostatic models the basic parameters are the classical effec- 
tive temperature T e g and gravity log(g). The effective temperature is run 
from 450kK (= 450,000K) to lOOOkK, in 50kK steps for the hydrostatic mod- 
els. In the grid of nova-type hybrid models the effective temperatures are 
550kK, 600kK and 650kK (extension to lower and higher values is work in 
progress). This range is based on comparison with observations from novae 
in a supersoft X-ray state (see section 5.4). 

The log(<7) is varied along with T c q, accounting for the fact that the 
radiation pressure increases with the luminosity and thus with T e g. In a 
hydrostatic atmosphere the effective gravitation determines the radial ex- 
tension (see chapter 3). Therefore, in order to compare atmospheres with a 
comparable effective gravitation, the following relation [Mih78] was used to 
define a default log(g) 



from which then small systematic deviations were made in both directions. 

In a later stage, it was found from the original derivation of relation (5.1) 
[Und49] that the right hand side is a factor of 10 too high. The right version 
is 



By that time a major part of the grid computations was already completed. 
From then on, only the lower log(g) parts of the grid have been continued, 



g > 65(T cff /10 4 ) 4 



(5.1) 



g > 6.5(T cfr /10 4 ) 4 



(5.2) 
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so that in the grid most models are computed around g being half of the 
(wrong) lower boundary given by equation (5.1). Consequently, all models 
shown in this section feature a rather high log(g). 

The results of these grids are shown in the figures 5.1 and 5.2 for the 
hydrostatic models and 5.3 to 5.8 for the hybrid models. The figures show 
the effects of the following parameters to the model spectra and temperature 



structures. 
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Figure 5.1: The effective temperature T e g is varied and a corresponding approximate 
effective gravitation is kept fixed using equation (5.1). The upper graph shows the spectra, 
the lower graph the temperature structures on the optical depth reference scale. 
A systematic change in the spectra and the temperature structures can be observed from 
changing the parameter T c g. 
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Figure 5.2: The log(g) is varied for two fixed effective temperatures. The temperature 
structures are plotted logarithmically against the density. 

There are two branches in the temperature structures in the tenuous part of the at- 
mosphere, corresponding to the two effective temperatures, and three branches can be 
discerned in the high density region, which correspond to the three values of log(g) that 
were used. Apparently, in these models the impact of these two parameters varies with 
depth in the atmosphere. In the 600kK spectra the log(g) affects the strength of the C VI 
ionization edge, which differentiates the spectra. In the 500kK spectra the influence of 
log(g) is very small. 
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Figure 5.3: In this figure the effective temperature T e g is varied and a corresponding 
approximate effective gravitation is kept fixed using equation (5.1). Also, all velocity field 
parameters are identical. The temperature structures are a bit more variable than in 
the hydrostatic models. Around a density of p = 10~ 9 g/cm? a dip develops with rising 
temperature. But there is a gradual development in the spectra with T e g, comparable to 
that for the hydrostatic spectra of figure 5.1. 
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Figure 5.4: The log(g) is varied for two fixed effective temperatures and velocity field. 
The temperature structures bundle in two in the tenuous part of the atmosphere, corre- 
sponding to the two effective temperatures, and in two other branches in the high density 
region, which correspond to the values of log(g). The same effect was observed for the 
hydrostatic models (figure 5.2). The differences in spectra caused by log(g) are the largest 
between 14.2 and 25. 3A (due to the Oviii and Cvi edges respectively), which also corre- 
sponds to the hydrostatic case. 
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Teff = 600kK, lg = 8.6, beta = 1.0, v_inf = 2400km/s 
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Figure 5.5: The mass loss rate M is varied for a fixed Voo and f3. This parameter 
obviously has a large influence on the spectrum, increasing the short wavelength flux. 
A gradual development can be observed longwards of the 14. 2A O viii ionization edge, 
in the continuum and in the lines. Also, the temperature structures show a systematic 
behavior. With increasing mass loss rates the density of the wind and thus the optical 
depth at the bottom wind layer (marked with a large cross) increases. Outwards of 
r rc f = 10 -5 no systematic temperature development can be expected (see section 4.8.1). 
The models are stable for the whole range of mass loss rates tested. 
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Figure 5.6: The /3 velocity field parameter is varied for a fixed mass loss rate and 
terminal velocity. This parameter has a large influence on the spectra. 
With a larger ft the velocity increases more slowly outwards, and from the continuity 
equation (3.3) it follows that then also the density falls more slowly. Therefore, the 
optical thickness of the expanding shell is larger. This can be seen in the lower graph, 
where the optical depth of the bottom layer of the wind (marked with a large cross) is 
larger for the higher values of /3. With the larger optical depth also the temperature at the 
transition layer increases through a stronger backwarming effect, and thus the spectrum 
appears hotter. 
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Figure 5.7: In this figure the asymptotic velocity Vca is varied and the mass loss rate M 

and P are kept fixed. This seems to have major impact on the spectra: higher velocities 

result in decreased flux in the short wavelength range, so the spectrum appears cooler. In 

the lower graph it can be seen, that indeed the temperature in the optically important 

region 10 -3 < r rc f < 10° is lower for the higher wind velocity models. 

This behavior can be understood in the same way as the plots in figure 5.6. An increase in 

the velocity for a fixed mass loss rate leads to an increase in the density according to the 

continuity equation (3.3). The opacity scales (roughly) with the density. As a consequence 

of the increased opacity then the temperature rises (backwarming effect). 

Compare also figure 5.8, where M is scaled linearly with v^o, so that the density structures 

are constant. 
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Figure 5.8: In this figure the terminal velocity Voo is changed, but at the same time the 
mass loss rate M is scaled by the same amount, so that the resulting density structures 
for the expanding envelopes are identical (see equation (3.3)). 

Compared with the differences due to only Voo variation (see figure 5.7) the differences 
between the models in this figure are very small. The spectra are very similar, see figure 
5.9 for a zoom in on a smaller part of the spectrum. Also, the temperature structures are 
very similar, except for the outer region that is optically very thin with r re f < 3 • 10 -5 . In 
those regions the wavelength averaged opacities become dominated by the contributions of 
single strong atomic lines, and therefore also the radiative equilibrium and the temperature 
are. This effect is called line cooling /heating, see section 4.8.5. The velocity field does 
influence the line cooling/heating effect, and therefore also the temperature structure. 
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Figure 5.9: This figures shows the same spectra as in figure 5.8, but now zoomed in to 
a smaller wavelength range and on a linear wavelength scale. 

In the spectra the continuum fluxes are almost the same. There is a stronger blue shift in 
the absorption lines for the models with higher velocities. Also, the details (narrow lines) 
seen in the low velocity spectrum (black curve) are smeared out with higher velocities. 
Therefore, the spectrum for the highest velocity in this example (green curve) is very 
smooth. 



5.2 LTE vs. NLTE 

The major part of the theory and methods described in sections 2 and 4 
dealt with the problem of computing NLTE models. The extra complica- 
tions in atmosphere modeling (in both methodology and computation time), 
introduced by relaxing the assumption of LTE, are yet to be justified. 

Whether LTE is a useful approximation for a specific atmosphere or not 
can only be ascertained from the direct comparison with a NLTE model. 
If the deviations from NLTE are small, then LTE is a good approximation, 
otherwise it is not. This direct comparison is done here for the expanding 
nova type and hydrostatic type models. Since the nova type models are 
based upon a hydrostatic core, a pure hydrostatic atmosphere is recovered 
in the limit for tenuous winds, i.e. very small mass loss rates. 

5.2.1 Expanding models 

In figure 5.10 a LTE spectrum is shown in comparison with a NLTE spectrum 
for a typical radially extended nova atmosphere of the type described in 
chapter 3. If a LTE spectrum is computed from the NLTE atmospheric 
structure, so that the atmospheric conditions for the models are identical 
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Figure 5.10: Here the spectra are shown for a typical nova-type model for four different 
situations. The models were computed in descending order (as shown in the legend), 
starting from the 'Pure NLTE' model (black), each of them using the previous model 
result as initial condition. The first two spectra (black and red) are computed from the 
NLTE atmospheric structure. For these two all atmospheric conditions were identical, 
except for the population numbers, being in statistical equilibrium (black) or in LTE 
(red). In the green model the temperature was released to adapt to the LTE opacities, 
but the density structure was kept fixed. The blue model is the pure LTE case, where the 
whole structure is released to adapt to the LTE opacities, i.e. the temperature, density 
and all derived quantities (like partial pressures). 

For the radially extended nova atmospheres treated in this work obviously LTE is a poor 
approximation. 



except for the population numbers, then the result is already very different 
from the pure NLTE spectrum. If only the temperature structure is released 
to adapt to the LTE opacities, then the resulting spectrum becomes already 
more different. If the atmospheric structure is released to adapt to the LTE 
opacities instead of using the NLTE structure, then the resulting spectra 
become even more different, see figure 5.10. The temperature structures for 
these two cases are plotted in figure 5.11. Obviously, the NLTE treatment 
is essential for this type of models. 
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Figure 5.11: These are the temperature structures for the different NLTE and the LTE 
spectra compared in figure 5.10. The NLTE temperature structure for the black and red 
spectrum coincide. With the LTE structure (blue curve) the outer parts of the atmosphere 
are already optically thick at very low densities, because the low ionization stages these 
temperatures and LTE yield very effective X-ray absorbers. 
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Figure 5.12: These spectra are for a hydrostatic atmosphere. The structures for these 
spectra are identical except for the population numbers. The most obvious difference 
is the strength of the ionization edges, a result that was also found for the nova-type 
atmospheres in figure 5.10. This is mainly caused by the ionization balance shown in 
figure 5.13. 



5.2.2 Hydrostatic models 

Whereas the departures from LTE are large for the hybrid nova-type models, 
this is not necessarily the case for other types of models, like the hydrostatic 
ones examined here. These hydrostatic atmospheres are very compact in 
contrast to the radial extension of nova-type models. In figure 5.12 a LTE 
spectrum is shown in comparison with a NLTE spectrum. The differences 
between LTE and NLTE are much smaller than for the nova-type atmo- 
sphere in figure 5.10. However, the discrepancies are still large enough to 
disqualify LTE as a good approximation, the more since the underlying 
structure was computed in full NLTE. So even if a good structure is avail- 
able, then LTE is still a poor assumption. The most striking difference in the 
spectra is the strength of the C VI ionization edge at 25. 3A. The ionization 
balances for LTE and NLTE for this hydrostatic model are shown in figure 
5.13. Clearly, C vi is suppressed by NLTE effects in the outer regions of the 
atmosphere. But differences are found not only in the ionization edges but 
also the detailed line features. This becomes visible when zooming in to a 
smaller wavelength range, which is shown in figure 5.14. 
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Figure 5.14: A zoom into a smaller wavelength range is plotted for the spectra of figure 
5.10. This displays that differences between LTE and NLTE for compact hydrostatic 
atmospheres are not limited to strong ionization edges but are also apparent in the line 
features. Note that the LTE spectrum is computed from the NLTE structure. From pure 
LTE much larger differences can be expected (see also figure 5.10, where this was shown 
for nova- type models). 
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5.3 Transformation of counts to flux spectra 

Where in the previous sections it was shown that systematic results are ob- 
tained from the new framework, the question if these models give a realistic 
description of nova atmospheres is yet open. In order to compare theoretical 
spectra with X-ray observations the interstellar extinction and instrumental 
characteristics must be considered. 

5.3.1 Interstellar extinction 

Accurate estimates of interstellar (IS) extinction are necessary for the anal- 
ysis and interpretation of almost all astronomical soft X-ray observations. 
The IS extinction in the X-ray wavelength range is far from gray. The soft 
X-rays are more strongly absorbed than the hard X-rays. Therefore, the 
IS extinction does not only influence the estimated bolometric luminosity 
(from a given radius and distance) of the object, but also the hardness ratio, 
i.e. the ratio of hard and soft photon counts. 

In this work, an IS model is used that is based on the absorption cross sec- 
tions compiled by [BM92] for 17 astrophysically important species. Specif- 
ically, the effective extinction curve is calculated with Balucinska-Church 
&: McCammon's FORTRAN code 1 , assuming solar abundances. The total ex- 
tinction is then proportional to the hydrogen column density. The model 
includes (continuum) absorption from bound-free transitions only. Although 
this code has not been updated since 1999, it is still the most widely used 
in literature. 

Also, a newer code exists, called tbabs [WAM00], that includes more 
fancy absorption features, like absorption lines from neutral elements. An 
example given by the authors that compares the absorption coefficients from 
an old and a new version of tbabs is shown in figure 5.15. Inclusion of these 
features is important when a chi-square method is used to determine the 
"goodness" of a fit. But the chi-square method poses some fundamental 
problems when used for IS absorbed spectra (see section 5.3.3), so that 
it is not used in this work and the fits are done "by eye". Apart from 
improving the fit goodness value x 2 -, the line absorption features provide 
a better basis for tuning the assumed interstellar chemical composition set 
from the comparison of observations with synthetic spectra, tbabs has not 
yet been used in this work. It is designed for XSPEC, a large spectral 
analysis package, and is not out-of-the-box suitable include in this work's 
code. Private communication with the authors and a detailed comparison 
of the IS absorption cross sections and its influences on the model fits is 
planned for the future. 

Whereas the slope of the 'Wien' tail of the spectrum The IS extinction 
inclines so strongly to longer wavelengths, that the slope 



1 Available at ftp://adc.astro.umd.edU/pub/adc/arcliives/catalogs/6/6062A/ 
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Figure 5.15: This plot compares the Ol ionization edge from a new version (red) of the 
interstellar absorption package tbabs with an oder version (black). This code is newer 
than commonly used code of [BM92], that also used in this work. Detailed comparison 
and possibly inclusion of this newer package into the code used for this work is planned 
for future. From: http://astro.uni-tuebingen.de/~wilms/research/tbabs/ 



5.3.2 Instrumental characteristics 

The 'high-res' X-ray observations used in this work come from the X-ray 
grating satellites Chandra and XMM-Newton, with the LETGS and RGS 
spectographs respectively. After the data reduction process, that deals with 
most of the instrumental characteristics, one obtains the number of photons 
that were detected per energy bin in the total exposure time. This reduction 
process was done by [Nes09] for all observations. The reduced data still 
contain higher order contributions or detector gaps that must be dealt with 
(see below). 

Given the effective area a of each bin i with bin width w and the total 
exposure time r, the photon count n can be converted to the observed 
(absorbed) flux F in the center wavelength A of the bin according to 

Ft = £ (5.3) 



with n bg being the background counts. 
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With the transmission coefficient T, being a function of the hydrogen 
column mass and the chemical composition xj, the unabsorbed flux becomes 

pahs 

- A (5.4) 



T x (N K , Xj ) 

The inverse transformation, from unabsorbed flux to counts, thus becomes 

ni = F Xi T Xi a Xi TW Xi ^- c + nY (5.5) 

The resolution of the synthetic spectra is significantly higher than that of 
the observed spectra. Therefore, the conversion process is done on the ob- 
served wavelength grid, and the synthetic flux points must be spread over 
that grid. For the comparison with observations, synthetic spectra must be 
transformed from the Lagrangian to the Eulerian frame, and for the Lorentz 
transformation no wavelength parallelization can be done. In order to save 
unnecessary computation time, the synthetic spectra are computed with not 
many more wavelength points than necessary, i.e. a spectral resolution 2 of 
roughly 3 times the observed resolution. This means, that in the spread- 
ing of the synthetic flux over the observed wavelength grid, not only the 
flux points that fall within the detector bin must be treated, but also the 
proportional parts of the next closest flux points on either side of the bin, 
according to their distance to the bin boundary. 



Chandra specific 

However, the Chandra LETGS spectra contain the overlapping contribu- 
tions of the first and higher diffraction orders. The relative strengths of 
the different orders are plotted in figure 5.16. The photons in the higher 
orders are measured in the spectrograph at the double, triple, etc. of the 
actual wavelength. Writing the photon counts in bin i as the sum of the 
contributions of the first three orders 

m - n^ g = n\ + n 2 + n\ (5.6) 

and using equation (5.5) one obtains 



Fx, 



he 



A* T x ,a x ,w x . 



tp rp 2nd A/2 

rii - F Xi/2 T Xi/2 a Xi w Xi/2 -^ 



3 rc j A/3 



(5.7) 



2 Note that the synthetic resolution can vary with wavelength. 
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Figure 5.16: This plot shows Chandra's higher order diffraction efficiencies relative to 
the 1st order efficiency. The dotted curves show the pre-launch model, the solid curves 
the current model predictions. From the Chandra X-ray Center homepage. 



This equation involves simultaneous values of F for different wavelengths. 
But due to the systematic wavelength shift in the orders, the equation can 
be solved recursively, starting at the small wavelength boundary. 

Disentangling the higher orders is very helpful in determining the proper 
IS extinction, parameterized by A/h, as the influence of Nn is the largest for 
the long wavelengths in which also the contributions of higher orders are 
most significant. It is stressed at this point, that the fits are very sensitive 
to A^h, so that a small change in A^h results in a large change in other fit 
parameters, like T e g. 

An example of the result from extracting the first, second and third order 
is shown in figure 5.17. The missing treatment of higher becomes apparent 
in the rise of first order flux for long wavelengths A > lOOA. 



V4743 Sgr, Mar 2003 (day 180) LETGS 7VH=1.3e21 cm" 2 
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Figure 5.17: The first 3 orders extracted from the spectrum of figure 5.f8. a\, a x n and af d are normalized, thus the relative strengths are 
unrealistic. 
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XMM-Newton specific 

The RGS spectra do not feature higher orders, but they have detector gaps 
(due to hot pixels and defective CCD chips), i.e. wavelength ranges in which 
the effective area is 0. These appear as strong spikes in the observed data. 
With a = the conversions between the representations cannot be applied. 
When the observed counts from the two independent detector channels, 
RGS1 and RGS2, are combined, then the number of gaps is strongly reduced 
and the quality of the spectra improves significantly. In order to cope with 
these gaps, both the effective area's and the observed counts are interpolated 
linearly. Apart from the gaps, still jumps in the effective area as large as a 
factor of 2 occur. 

In order to ease the interpretation the linear counts/bin spectra, these 
are converted to a constant effective area of a = 1, thus removing the pseudo 
absorption features in the observational data. That means that the counts 
are plotted as they would have been seen with an instrument that is equally 
sensitive over the whole wavelength range. 



5.3.3 The fitting procedure 

In order to compare the theoretical and observational data, the data must 
be transformed into one of both representations: unabsorbed flux or IS ab- 
sorbed counts/bin. In literature, the counts/bin representation is the most 
common. The advantage of this representation is, that the transformations 
are easiest to do in forward direction, like simulating (instead of disentan- 
gling) the overlapping higher orders of diffraction. But it is stressed here, 
that it is important to do the comparison between synthetic and observed 
spectrum in both representation simultaneously. 

The counts/bin spectra plotted on a linear scale reveal the statistical 
nature of low counts detections, that are important in X-rays. It gives an 
immediate impression in which wavelength ranges the detection was strong 
or weak with respect to the background signal. Also, it allows to show the 
number of counts that are interpreted as background in the same plot. 

On the other hand, the unabsorbed flux, plotted logarithmically, is very 
sensitive to the IS extinction. As in the red tail of the (absorbed counts/bin) 
observation the actual flux from the object is attenuated typically by a factor 
1000, the signal that is left in the observer representation is very small, so 
that deviations of a factor of 2 or more do not catch the eye. In the flux 
representation, where the observed signal is corrected for IS extinction by 
the same large factor, even small deviations in the absorbed model become 
apparent. 
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'By eye' instead of chi-square 

Often it is argued, defending the applicability of the chi-square method as 
fitting procedure, that the reproduction by the synthetic spectra in the fit 
is the most important in those regions where the signal is the strongest. 
This is true for pure statistical significance tests, but not for comparing two 
physical models: while the reproduction of the observation is important 
in those regions where the original, unabsorbed flux would have been the 
strongest (the IS extinction has no physical relation to the actual flux of the 
object we want to model), also the reproduction of the overall shape of the 
spectrum, including the weak fluxes, is required for a physically acceptable 
model. 

Using the chi-square method has the very welcome effect, that Ah can 
be treated as an additional free parameter in the modeling process. This 
yields supposedly better fits in the absorbed counts/bin representation. A 
bad Nu in such a fit manifests in the red (absorbed) tail of the spectrum, 
where the slope of the model does not meet the slope of the observation. 
In this tail the model counts can be off by a relatively large factor, a factor 
that would certainly catch the eye in the region of maximum counts. On the 
other hand, constraining the Ah accurately from both the counts/bin and 
logarithmic flux plots, poses a strong restriction on the determination of a 
good fitting model spectrum, and leads to the disqualification of fits that 
'look better' (at first sight) in a linear counts/bin plot only with a wrong 
value of Ah- 
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5.4 First 'fits' to X-ray nova observations 

In sections 5.1 and 5.2 it was shown that with the improved theory and the 
new methods developed in this work (chapters 2 and 4) systematic results 
are achieved from the NLTE nova models. This is the prerequisite for useful 
modeling of observational data. 

Once this is made sure, the model spectra can be compared to observa- 
tional data, as the final goal of this work is to improve our understanding 
of the nature of the observed objects: the classical novae. The procedure 
to attain this goal is to construct a model, physically as realistic as possi- 
ble, and vary the assumptions in the model and compare the results with 
the observations. The assumptions in the model are being constrained by 
the generally accepted theoretical and observational results. The basic as- 
sumptions about the structure for nova models in this work are described in 
chapter 3. 

Another important assumption concerns the abundances of the chemi- 
cal elements in the model atmosphere. The abundances do not only have 
a major impact on the opacities, and thus on all model results, but also 
play an important role in the evolution of the nova. Sophisticated nova 
evolution models exist that can simulate nova outbursts, the TNR scenario 
(see section 1.1). These depend sensitively on the abundances of the pre- 
outburst configuration. With such models the nucleosynthetic yields can 
be traced throughout the outburst, and thus detailed chemical composition 
sets can be derived for the ejected material, see [SST74], [PSS78], [PSS79] 
and [SIH+09]. 

However, in the SSS phase the ejecta become optically thin and deeper 
layers become visible. These have a history that can be very different than 
the history of the matter ejected during the initial outburst. In general, there 
are three different sources of material with different chemical compositions. 

1. Material that is accreted from the companion, which is commonly 
assumed to have solar abundances. 

2. Material that is dredged up from the white dwarf. This material can 
be either CO enhanced, or ONeMg enhanced. 

3. Previously thermo- nuclear ly processed material from the nova that is 
reaccreted after the outburst. This material is strongly NO enhanced. 
The C abundance of the ejecta depends on the white dwarf material, 
but even novae with a ONeMg white dwarf are expected to have a 
super solar C abundance. 

Being a mixture of these components, the chemical composition of the deep 
layers that become visible in the SSS phase is not unlikely to be different 
from the composition sets available in literature, which are derived from 
either theory or observations for earlier stages after the outburst. 
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The SSS phase allows a 'deep' insight in the inner regions of the nova at- 
mosphere, 'close' to the region where nuclear burning takes place. Therefore, 
it will be very interesting to accurately determine the chemical compositions 
that are 'seen' in this phase. 

However, variation of the chemical abundances introduces a large number 
of additional free parameters to the model. And given the computational 
requirements per model, the spectral analysis (which involves fine-tuning the 
abundances) needs to be done in a systematic step-by-step way. The 'bulk 
grid computation' method used in section 5.1 to explore the parameter space 
of the basic atmosphere structure parameters, is practically not suitable 
because the dimensions of such grids become too large. Detailed spectral 
analysis has not yet been done in the scope of this work but is planned for 
the near future. 

The models that are compared to the observations in this section all 
use solar abundances [AGS + 05]. It is a standard set of abundances, that is 
suitable to test the models, but is possibly not very realistic. 

5.4.1 Expanding models 

In figures 5.18 to 5.27 the synthetic spectra for the expanding nova type 
models computed so far are compared with all well exposed 'high-resolution' 
observations of novae in a supersoft X-ray state available to date: 
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V4743 Sgr, February 2004, day 526, LETGS 
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RS Oph, March 2006, day 39.7, LETGS 
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RS Oph, April 2006, day 54.3, RGS 
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RS Oph, April 2006, day 66.9, LETGS 
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V2491 Cyg, April 2008, day 39.9, RGS 
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.27: 


V2491 Cyg, April 2008, day 49.7, RGS 



The V4743 Sgr models have a terminal velocity of = 2400 km/s, the 
RS Oph models have Voo = 1200 km/s, and the models for V2491 Cyg have 
Vqc = 4800 km/s. These velocities have not been tuned. They originated in 
the grids made to explore the parameter space, shown in section 5.1, as the 
half and the double of 2400 km/s, which was originally chosen as default 
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value with the observed line shifts for V4743 Sgr [NSB+03] in mind. All 
models have /3 = 1.5. 

The distances to the objects are assumed to be 3.9, 1.5 and 10.2 kpc for 
V4743 Sgr, RS Oph and V2491 respectively. These distances only affect the 
fit parameter i?wD that scales the synthetic flux to the observed flux. 



V4743 Sgr, Mar 2003 (day 180) LETGS 7VH=1.3e21 cm" 2 
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Figure 5.18: Expanding model 




Figure 5.19: Expanding model 



V4743 Sgr, Jul 2003 (day 302) LETGS A/H = 1.3e21 cm" 2 




Figure 5.20: Expanding model 




Figure 5.21: Expanding model 




Figure 5.22: Expanding model 




Figure 5.23: Expanding model 




Figure 5.24: Expanding model 




Figure 5.25: Expanding model 



V2491 Cyg, Apr 2008 (day 39.9) RGS NH = 3.3e21 cm" 
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Figure 5.26: Expanding model 




Figure 5.27: Expanding model 
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5.4.2 Discussion: Expanding 

In the comparison between the models and the observations all X-ray grating 
(Chandra and XMM-Newton) spectra were used that are available to date. 
For none of the observations the models have been fine tuned, but rather it 
has been shown, that for all available observations a model can be found that 
at least reproduces the general shape and some detail of the observation. 

For the interpretation of the fit quality, it must be kept in mind, that 
the steps in the values of the basic model parameters in the grid were large. 
In the process of finding the 'best fit' model, this often results in one model 
being off in one direction, and the next model in the grid being off in the 
other, so that an increase in the parameter resolution in the grid would 
certainly improve the quality of the 'best fit'. 

V4743 Sgr 

The V4743 Sgr models have a terminal velocity of Vqq = 2400km/s. In 
all of them the flux in the blue tail of the model spectra is too strong. 
This is most apparent shortwards of 22. 5A, which exactly matches the Nvi 
ionization edge, so that increasing the N abundance to super solar could fix 
this problem. The models with T c g = 500kK for the September and April 
observations show a somewhat too strong emission component in some of the 
P Cycni profiles (around observed wavelengths of 27. 4A and 30. 2A). On the 
other hand, for other lines these emission components fit the observations 
very well (around observed wavelengths of 26. 8A and 32. 6A) 

RS Oph 

In all RS Oph models, that have a terminal velocity of = 1200km/s, 
but especially in comparison with the last two (April) observations, some 
absorption lines are too narrow and too deep. This is probably due to 
the neglected instrumental response matrix. The instrumental accuracy 
is usually not exactly as high as the resolution, as photons of a certain 
wavelength are smeared out over a few bins in the detector. Such a smearing 
effect would make the lines in the model a bit broader and less deep. 

The middle observation, of day 54.3, is much smoother than the first and 
the third observation. The spectral features in the model are also found in 
the observations, but there they are much smoother. For this observation 
a periodic variability was discovered [OPG + 06]. It would be interesting to 
divide this observations in multiple time frames, in order to see if there is a 
temporal development in the observation that causes details to be smeared 
out in the time integrated data. 

For the day 54.3 and day 66.9 observations, when going from lower to 
higher wavelengths there is an obvious step up (increase) in the observed 
flux at 18.6A, being the ionization edge of Nvn. This is not followed by 



124 



CHAPTER 5. FIRST RESULTS 



the models. Again, increasing the N abundance would probably fix this 
problem. In the first model, for day 39.7, this problem doesn't seem to 
exist. It is hard to believe, that the abundances would have changed very 
strongly since then. It will be interesting to find out if an abundance set 
that suits the later observations will also supply models that fit the March 
observation. 

V2491 Cyg 

V2491 Cyg is the fast expanding nova. The models were computed with 
Voo = 4800km/s, twice the speed as used from V4743 Sgr. For this outburst 
only two grating observations exist. And these observations show very large 
differences. In contrast to the other two sources, where the bolometric lu- 
minosity (displayed in the legend of each plot) is almost constant for the 
observations, the luminosity decreases by a factor of three. The first obser- 
vation, on day 39.9, was done at almost maximum X-ray brightness (see the 
lightcurve in figure 1.2). After maximum brightness, the count rate quickly 
declined. 

In the model for the first observation, with T c g = 700kK, the Ovu 
ionization edge at 16. 8A is too weak. Also, the characteristic O absorption 
lines at 19.0A and 21. 6A are much too weak in the model. However, in the 
cooler model for the second observation the edge and the lines are strong 
enough. The influence of changes to the abundance can be quite different for 
such different models, but the difference needed here is very large. Probably 
not only the abundances need to be tuned but also the model structure. 

Furthermore, for both observations the N vi edges at 22. 5 A are too weak 
and also the Nvn line at 24. 8A and the Nvi line at 28. 8A. 

The Cvi absorption edge at 25. 3 A is too strong in the model for day 
39.9 and a bit too strong in the model for day 49.7. 

5.4.3 Hydrostatic models 

In figures 5.28 to 5.37 the spectra from hydrostatic models computed so far 
are compared with all well exposed 'high-resolution' observations of novae 
in a supersoft X-ray state available to date: 

• figure 5.28: V4743 Sgr, March 2003, day 180, LETGS 

• figure 5.29: V4743 Sgr, April 2003, day 196, RGS 

• figure 5.30: V4743 Sgr, July 2003, day 302, LETGS 

• figure 5.31: V4743 Sgr, September 2003, day 371, LETGS 

• figure 5.31: V4743 Sgr, February 2004, day 526, LETGS 



• figure 5.33: RS Oph, March 2006, day 39.7, LETGS 
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• figure 5.34: RS Oph, April 2006, day 54.3, RGS 

• figure 5.35: RS Oph, April 2006, day 66.9, LETGS 

• figure 5.36: V2491 Cyg, April 2008, day 39.9, RGS 

• figure 5.37: V2491 Cyg, April 2008, day 49.7, RGS 

Since the models do not treat the expansion, the spectra are blue shifted 
with 1000, 2400 and 4000km/s for the V4743 Sgr, RS Oph, and V2491 Cyg 
observations respectively. 



V4743 Sgr, Mar 2003 (day 180) LETGS 7VH=1.4e21 cm- 




Figure 5.28: Hydrostatic model 




Figure 5.29: Hydrostatic model 




Figure 5.30: Hydrostatic model 




Figure 5.31: Hydrostatic model 




Figure 5.32: Hydrostatic model 




Figure 5.33: Hydrostatic model 
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Figure 5.34: Hydrostatic model 




Figure 5.36: Hydrostatic model 




Figure 5.37: Hydrostatic model 
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5.4.4 Discussion: Hydrostatic 

Just like for the expanding models shown in section 5.4.1, the hydrostatic 
models are compared to all X-ray grating (Chandra and XMM-Newton) 
data that are available to date. The models have not been fine tuned, they 
all use the same solar abundances as the expanding models. The parameter 
resolution in the grid (e.g. the step size in T e g- is 50kK) is identical with the 
expanding models grid. 

V4743 Sgr 

The hydrostatic models reproduce the observations of V4743 Sgr quite well, 
in this respect the quality is comparable to the expanding models, except 
for the last observation, day 526. The model with T e g- = 450kK features a 
strong Cvi ionization edge at 25. 3A that is not present in the observations 
and very weak in the expanding models. 

However, the reproduction of the detailed spectral features is much worse 
than the expanding models. There is hardly any absorption feature that 
matches the observation in line shape, line strength and line center. The 
ad-hoc blue shift does move the lines closer to the observed line centers but 
since lines are formed in different regions in the wind with different local 
velocities, the global constant shift is not accurate. 

RS Oph 

The hydrostatic models compared with RS Oph all show a much so strong 
O vii ionization edge at 16.8A, so that shortwards of this edge the flux much 
too weak. Surely, this could be compensated with a reduced O abundance, 
but it is interesting that a sub-solar O abundance is not required for the 
expanding models for RS Oph. Also, for all observations the O VII absorption 
line at 21. 6A is much too strong in the model spectra, which supports the 
idea that the O abundance should be decreased in order to improve the fit. 

V2491 Cyg 

Just like for RS Oph, the hydrostatic models show a much too strong O vii 
ionization edge and too strong O absorption lines, like the Ovu line at 
21. 6A. Decreasing the O abundance by a factor of 10-100 would certainly 
improve the fit. But again, this is not found for the expanding models. 
On the contrary, the O abundance was concluded to be too low in for the 
expanding models. 

Furthermore, for both observations the N vi edges at 22. 5 A are a bit too 
weak and also the Nvn line at 24. 8A and the Nvi line at 28. 8A. This was 
also found for the expanding models. 
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The Cvi absorption edge at 25. 3 A is too strong in the model for day 
49.7 and a bit too strong in the model for day 39.9. Also, in the expanding 
models a too strong C vi edge was found, however there the effect was clearly 
stronger for day 39.9 instead of for day 49.7. 
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Chapter 6 

Conclusion and Outlook 



In this work nova atmospheres are treated as a hybrid type atmosphere 
consisting of a hydrostatic core with an expanding envelope on top of it, 
see chapter 3. In section 4.2 it was shown that these nova models need to 
be treated as pure NLTE models. This requirement introduces a number 
of complications that have lead to new methods and improvements to the 
code. In section 5.1 it was shown that systematic results are obtained from 
the new framework. This is a prerequisite for doing any useful modeling 
of observed spectra. However, this is not yet a guarantee for successful 
modeling of observations. Therefore, the next step has been to compare the 
model results with observations in section 5.4. 

It is found, that for all well exposed grating observations available to 
date models can be found in the test model grid that match the observations 
quite well, section 5.4.1. Also, comparison of the observations was done with 
hydrostatic models that were computed with the new version of code, section 
5.4.3. Generally, the general shape of the observations can be reproduced 
reasonably well, but they are by far inferior to the expanding models in 
matching the detailed spectral features. A very interesting finding from the 
two sets of comparisons is that often contradicting conclusions are drawn 
from the expanding and hydrostatic models about the corrections that need 
to be made to the abundances in order to improve the fit, sections 5.4.2 and 
5.4.4. 

The new framework appears to be ready and working. The first results 
from the new models show remarkable agreement with the observations. 
With the final goal to understand the physical nature of novae, this is a 
promising start in the important task of determination of the physical con- 
ditions of nova atmospheres in the late SSS stage. But a lot of work is yet 
to be done and numerous improvements and refinements to the models are 
outstanding. 
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6.0.5 Observational objectives 
Tuning the models 

• The parameter resolution (stepsize) used in the model grids is yet 
poor, so when the observation falls right between two adjacent models 
then the fit is far from optimal. Therefore, the grid resolution must 
be increased. 

• All models computed thus far used solar abundances. Although the 
precise chemical composition that is 'seen' in this late SSS phase of 
the nova outburst is not yet well determined, it can be expected from 
measurements of earlier stages in the outburst that solar abundances 
are not very realistic. In the process of fine tuning the abundances pre- 
sumably also the values of the basic atmosphere structure parameters 
will change, so that this is an iterative process. 

• In this work no attention has been payed to the details of the atomic 
data. As shown in previous work on X-ray nova with PHOENIX [Pet05] 
the atomic data can have a major influence on the spectrum. However, 
the new models of this work significantly differ from the models in the 
previous work. Therefore, the conclusions drawn must be tested and 
revised where needed. 

• The interstellar extinction model used in this work is rather basic. 
Newer, more fancy IS models exist in literature [WAM00], see also 
section 5.3.1. The influence of the newer model must be analyzed. 
Given the very large impact of the IS extinction, improvements in 
this field could have an effect on the fit quality of synthetic spectra to 
observations. And, more importantly, in the process of tuning a model 
to an observation errors from a wrong IS Model could lead to wrong 
fit parameters. 

• Another important thing to do is to determine the uniqueness of a 
good fit, or in other words, to determine the 'error bars' on the fit 
parameters. 

Interrelation with nova evolution models The results of good model 
fits could provide useful constraints for the development of hydrodynamical 
nova evolution models, like [SIH+09]. And those, in turn, could provide 
constraints for the radiative transport models. Collaboration in this field of 
work looks very prospective. 

Fragmentation of long exposure times Many of the grating spectra 
are exposed well enough to allow for segmentation of the data into smaller 
time frames. Each time frame can then be analyzed separately, as done for 



141 



example in [NSB + 07]. Several observations show puzzling short timescale 
variations in the lightcurve. Investigating these variations would be inter- 
esting to begin with. 

New observations Although the grating satellites are already coming 
to age, it is not easy to get observation time for X-ray novae. The new 
models presented in this work will provide the means to analyze and interpret 
the observations. This is a very important argument in the race against 
competing research field, and due to the lack of this ability it has become 
very hard to get allocation time for nova observation. 

Other observational data Apart from grating spectra, there is a huge 
amount of low-resolution X-ray spectra from the Swift satellite. Although 
these data do not provide much the spectral detail, they could still be valu- 
able, especially when grating spectra are available for the same object. The 
grating spectra could be used to determine a good fine tuned model and 
the Swift spectra could provide the information to determine changes in the 
atmosphere from that point. 

Multi-wavelength band studies Since X-ray observations involve some 
typical difficulties, e.g. the uncertainties introduced by unconstrained IS 
extinction, extension to multi-wavelength band studies could provide very 
interesting new insights. Swift does do simultaneous multi-wavelength band 
observations already. 

6.0.6 Theoretical objectives 

Apart from the numerous objectives that involve comparison with observa- 
tions, there is a number of outstanding theoretical questions that need to 
be elaborated further. 

• The line profiles used in the models must be sophisticated. As de- 
scribed in section 2.5.2 the stark broadening treatment as currently 
implemented is not accurate for the atmospheric conditions close to 
the white dwarf boundary or below. 

• In section 4.5 a method has been proposed that combines the linear and 
parabolic interpolation of the source function over the characteristic 
rays, between the intersection points with the concentric model layers. 
As shown, that method will improve the accuracy of the solution of 
the radiation field. This has two advantages. In the first place, this 
could save idle iterations because the consistency between radiation 
field and radiative rates will be reached more quickly. Secondly, this 
method should reduce the number of ALI iterations required to reach 
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the prescribed accuracy. On itself the radiation transport part of the 
code is computationally very quick, but this part of the code must 
be computed sequentally for all wavelength points. Therefore, when 
scaling the code to many processors, the radiative transport is in the 
end the bottleneck. Speeding up the small radiative transport piece 
of code with the new interpolation method would yield a large overall 
performance gain. 

• The target flux for static atmospheres is constant with radius. For 
moving atmospheres the target flux can be transformed to the La- 
grangian frame with the equations given in section 2.11.2. However, 
the gas also cools when it expands. This effect has silently been ig- 
nored in the derivation of the approximate non-static generalization 
of the UL temperature correction method. Also, the dJ^o/dr term in 
(2.119) that vanishes in the static case has not been discussed. These 
effects could have a significant impact on the temperature structure, 
and thus on all model results. The impact depends on the velocity 
and density fields of the wind. 

• It was shown in section 4.8.1 that in NLTE the applicability of the UL 
temperature correction method is limited on physical grounds. The 
method is based on the assumption that the temperature of the gas 
can be derived from the radiation transport equation. This can only 
work if the opacities, that is where the properties of the gas enter, 
depend significantly on the temperature. In case of LTE that is given. 
In the other extreme of large NLTE effects it is not, and the UL method 
fails. 

In literature a completely different method is known that is based on 
the thermal balance of electrons [KPP99]. This method has already 
been implemented in the code, but could not yet be tested or used due 
to time constraints. For the nova models the temperature in the outer 
regions is does not have a large influence on the models (see section 
4.8.4 for a discussion), but (also for other types of models) this might 
yet be interesting to test and compare with the current temperature 
methods. 

6.0.7 3D models 

The ID spherical symmetry that has been assumed in this work clearly 
is a poor approximation given the binarity of a nova system. However, 
at present there are no 3D hydrodynamical models for novae that could 
provide a realistic atmospheric structure for 3D radiation transport. But 
the transition to 3D radiative transport models will be important as soon 
as such 3D structures become available. 
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6.0.8 Application of the new methods to supernova models 

Finally, the framework described in this work has been developed for X-ray 
novae. But none of the methods is limited in application to this specific class 
of objects. On the contrary, most of the improvements will help as much in 
other types of models. Especially supernova models will be very interesting 
to model given their similarities with nova models: high expansion velocities, 
large radial extension and strong NLTE effects. The more, since it has never 
been possible before to treat all species in pure NLTE, including all relevant 
ionization stages of Fe, Co and Ni (see sections 4.6.5 and 4.6.6). 
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Appendix A 

Problems with broad line 
theory 

Radiation can equivalently be described on the wavelength scale or on the 
frequency scale. The radiation field has a certain energy distribution that 
is interpreted in dA units or in d^ units. However, one has to be careful 
with assumptions that are made in one representation or the other. Espe- 
cially when dealing with broad atomic line profiles, a number of problems 
arise. Before discussing the bigger problem with the Einstein coefficients, 
line profile symmetries are discussed, as an example of the difference between 
assumptions in different representations. 



A.l Line profile asymmetry 

If a line profile is symmetric around the line center in the wavelength repre- 
sentation, it is not symmetric in the frequency representation, and vice versa. 
This is caused by the non-linear transformation between the two represen- 
tations. For very narrow line profiles, like Gaussian Doppler profiles, this 
asymmetry is a negligible effect. However, for broad profiles, like Lorentz or 
Voigt profiles (equation (2.64) or (2.57)), this asymmetry becomes apparent. 

The frequency-Lorentz profile, equation (2.64), converted to the wave- 
length representation is given by 

c A^/vr c 



x2 (f-i) 2 + (w A2 



(A.l) 

[Au d XX /c]^ 1 AX d 



it (A - A) 2 + [A^AAq/c] 2 vr (A - A) 2 + (AA d A/A ) 2 
where in the first step equation (2.32) was used and in the last step the 
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Lorentz width parameter converts as 

Au d = -^AX d (A.2) 
A o 

This profile is still area normalized but not symmetric around Ao in the 
wavelength scale. 

The wavelength-Lorentz profile is given by 

1 AX d 

^ = ^(A -A) 2 + (AA d)2 (A - 3) 

and comparison with the converted frequency-Lorentz profile, equation (A.l), 
shows that they differ by a factor of A/ Ao in one of the terms in the denom- 
inator. This factor is approximately 1 when close to the line center, but for 
the line wings these factors significantly influence the profile shape. 



A.2 The Einstein coefficients B 

A very important difference between the wavelength and the frequency rep- 
resentation lies in the definition of the Einstein coefficient for absorption and 
induced emission B. In short, for broad lines, if B is defined constant in the 
frequency representation, then it is wavelength dependent in the wavelength 
representation and vice versa. In section 2.5.2, where the bound-bound rates 
and opacities were developed, B was defined per intensity in the frequency 
representation I u , equation (2.27). In order to see the difference to B de- 
fined per intensity in the wavelength representation I\ here the expression 
for the upwards rates Rij are derived in the wavelength representation. 
The absorption rate per particle, in analogy to equation (2.27), is 



Rij — 



poo 

/ 0xJxd\ (A A) 

Jo 

where Bf- is the Einstein coefficient defined as the rate per units of inten- 



conversion rule for B^ follows 



sity in the wavelength representation. Using equation (2.33) the following 
fc 

Bf 3 = B^ (A.5) 

This means that if By is defined constant with A then Bf- is not, or if Bf- 
is defined constant then Bij is not. 

However, the Einstein relations derived in the wavelength representation, 
in analogy to equation (2.42), become 

Aji = ^jr^ (A.6) 
B% = gi/gjB* (A.7) 
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Note that Aji is independent of the representation in which the radiation 
field is described. Using equation (2.43) a different conversion rule for 
follows 

Bl 3 =B^ (A.8) 

The contradiction between equations (A. 5) and (A.8) results from the 
assumption of narrow lines made in the derivation of the Einstein relations. 
If lines are infinitely narrow, then the contradiction is resolved, as the rates 
and opacities are only evaluated at the point A = Ao- But for broad lines, ob- 
viously, the assumption of narrow lines made in the derivation of the Einstein 
relations is not consistent with the conversion of the frequency representa- 
tion result to the wavelength representation. Since different representations 
of the radiation field are physically equivalent, and in the conversion be- 
tween them no additional assumptions are made, the conclusion drawn here 
is that for broad lines the Einstein relations are not valid. 

In some literature, e.g. [Mih78], the Einstein relations are derived as 

= ^B JZ (A.9) 

Bji = gi/gjBij (A.10) 

However, this is inconsistent with using the Boltzmann level distribution to 
write the ratio of the TE level populations n« and rij in equation (2.42) as 

TE 

!h = 9i e hu lj/k T ^ 9i e hu/kT (A n) 

nj E 9j 9j 

The reason why the first Einstein relation is written in the form of equation 
(A.9) is to impose a uniformity in the expressions of opacities and rates 
between the bound-bound and bound-free cases. For the latter indeed pos- 
sesses the factor 



A. 3 Wavelength vs. frequency representation 

The physically right way to deal with this problem would be to separately 
derive expressions for the A and B coefficients, instead of relating them 
using equilibrium arguments. However, these would depend on the precise 
physical conditions of the system. Generally, they should have to be derived 
from a fully quantum mechanical description. But such treatment is already 
quite complicated for Hydrogen, see [Mih78], let alone for much larger atoms 
as oxygen or even Iron. This goes far beyond the scope of this work. And 
data of this quality are not available at the present time. 

But even if the Einstein relations are not exactly right, they might still 
give useful approximations. In this work, the two variants of equations (A. 5) 
and (A.8) were examined. 
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A. 3.1 Frequency variant 

The first variant was derived in the frequency representation in section 2.5.2. 
It uses frequency- Voigt profile functions, equation (2.57). Those are sym- 
metric around the line center on the frequency scale. The rates and opacities 
expressed in the A coefficients, that are independent of the choice of repre- 
sentation, for the frequency variant are 

f°° 4>\—J\ d\ (A.12) 
Ihc gi J Q c 

1+ 2tl < A - 13 » 

-Aji-rij^x (A. 14) 

o 



Rij 




Rji 


= Aji ( 
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Vij,x 


47rA ( 
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Xij,X 


47rA| 



A. 3. 2 Wavelength variant 

The second variant is based on the wavelength representation of the radia- 
tion field. The Einstein B coefficients are defined as the rates per units of 
intensity in the wavelength representation. In this variant wavelength- Voigt 
profile functions are assumed, which are symmetric around the line center 
on the wavelength scale. The expression for the wavelength- Voigt profile are 
similar to the frequency- Voigt. The dispersion width parameter converts as 
equation (A. 2), and similarly the Doppler width 

Av D = ^AA D (A.16) 



The u parameter becomes the dimensionless wavelength offset instead of the 
dimensionless frequency offset, but the Voigt damping parameter a, equation 
(2.61), is equivalent in both representations. 



i,(a,u)=i,(a,u) = ^^ (A.17) 

»=W (A - 18) 

a=£± = pL (A.19) 
AA D Az/ D 
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The rates and opacities expressed in the A coefficients, that are independent 
of the choice of representation, for the wavelength variant are 



A 5 „. ,00 



3 

2hc 2 Jo 



R% = An ( 1 + ^ / ip x J x dX ) (A.21) 



A he 



Vtj,x = ^x~ A n n 3^ ( A - 22 ) 

^=^SfX n ^- n ^) (A - 23) 



A. 3. 3 Flux comparison 

If the wavelength variant is used then a number of problems arise in the 
deeper atmosphere layers. The extinction of far line wings in the wave- 
length variant is almost independent of the wavelength, equation (A.23), 
since there the profile function is very flat, see figure A.l. The line profiles 
have wavelength symmetry, so they easily extend to small wavelengths 1 . 
At the same time, the energy of the radiation field per dA increases lin- 
early with decreasing A, so that the extincted energy from the far blue line 
wings increases linearly towards short wavelengths. In the deep layers of 
the atmosphere, the radiation field is approximately Planckian, and as the 
temperatures are high, the peak values are at short wavelengths. Because 
of the efficient far wing extinction in this wavelength range the bulk of the 
flux is strongly absorbed. Comparing the same layer in both variants shows 
that this effect decreases the integrated Eddington flux H in the wavelength 
variant by a factor of 4. 



A. 3. 4 Mean radiation field comparison 

Since the bound-bound contribution to the source function is only planckian 
in the line center, see section 2.5.2, the bound-free opacities are responsible 
for generating a planckian continuum source function. But in the wavelength 
variant the bound-free opacities do not dominate the bound-bound opaci- 
ties of the line wings at short wavelengths, because the blue line wings are 
so strong. Thus the source function in the continuum does not accurately 
approximate the Planck function, and so the mean intensity does not. In 
LTE the radiation field I\ is not exact planckian since then the Eddington 
flux H\ would be zero, like described in section 2.4. So in approximate LTE 

1 For a line with wavelength symmetry and the center at Ao = 100A extending the 
line from 10A to lA makes the wing just a factor of 0.1 broader, whereas for a line with 
frequency symmetry and the center at 1/0 = c/(100A) extending the line from c/lOA to 
c/lA makes the wing a factor of 10 broader. 
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Figure A.l: The top graph shows the source function S\, the related radiation moments 
J\ and H\ and the Planck function B\ for a deep layer in the atmosphere, layer 120 of 
128. In the two lower graphs the opacities are plotted, which generate the source function, 
equations (2.24) and (2.25). Here the wavelength variant of the bound-bound opacities, 
equations (A. 23) and (A. 22) are used. 

This variant causes trouble in the inner parts of the atmosphere. There the expected 
approximate Planckian shape of S\ is not reproduced. This is caused by the fact that the 
blue tails of the line wings are so strong that the bound-free opacities no longer dominate 
the continuum extinction. Furthermore, the bulk of the flux in deep layers is at short 
wavelengths and is therefore efficiently absorbed. Radiation transport shows, that this 
lowers the integrated Eddington flux H by a factor of 4 in comparison with the frequency 
variant. 

In figure A. 2 the same plots are shown for the frequency variant. 
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Figure A. 2: This figure shows the same quantities as figure A.l with the only difference 
that the bound-bound opacities are computed in the frequency variant, equations (A. 15) 
and (A. 14). 

This variant is superior to the wavelength variant. It produces an accurate Planckian 
source function above 3A, where the bulk of the radiative energy lies. Below 3A the wings 
of the bound-bound emission profiles are stronger than the bound-free emissivity, leading 
to a departure from the pure LTE source function. The only wavelength dependence in 
the emissivity comes from the profile functions. So if this departure of the source function 
is wrong, then the profile functions tpx are. The same effect occurs in the wavelength 
variant, of figure A.l. But since the radiation field is weak below 3A this does not have a 
significant influence on the models. 
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indeed deviations may exist. However, in the mean intensity J\ the devia- 
tions should be small, much smaller than those produced in the wavelength 
variant. In the frequency variant this problem does not occur, shown in 
figure A. 2. This is a strong indication that the wavelength variant of the 
expressions for the opacities, equations (A. 15) and (A. 14), are unusable. 

A. 3. 5 Departure coefficient comparison 

The inaccuracy in the mean intensity, together with the wavelength variant 
of the rate equations (A. 12) and (A. 13), leads to radiative rates that bring 
the population numbers out of their LTE balance in the innermost regions 
of the atmosphere. This is shown by the NLTE departure coefficients bi for 
all levels of oxygen in figure A. 3. With the wavelength variant, when going 
from the outer layer 1 inwards the departure coefficients converge towards 
bi = 1 up till around layer 100. From there on inwards they diverge again. 
With the frequency variant, the bi converge smoothly in the deepest regions 
of the atmosphere. 

A. 3. 6 Conclusion 

From the three problems that go along with the wavelength variant it is 
concluded that the description derived from the frequency representation 
and transferred to the wavelength scale is superior. This means that the 
Einstein coefficients B favorably are interpreted as constant in the frequency 
domain, so that they become wavelength dependent after conversion to the 
wavelength domain. 

A. 4 Frequency vs. original variant 

Originally in PHOENIX neither of the two variants described in the previous 
section were used. In the original variant used in PHOENIX the B^ are 
independent of wavelength and defined as the rates per unit of Intensity 
in the wavelength representation. Thus far, the approach is analog to the 
Wavelength variant of the previous section. 

But in the expressions for the opacities and rates the A ^ and B° { (where 
the o stands for 'original') are derived from the B\- using another form of 
Einstein relations 




(A.24) 





LTE 



" c he \ 

g A fcT g— XfcT B?) = 



(A.25) 



This form follows from equation (2.41) when the summations over all atomic 
energy states i and j of all atomic species are neglected. This means that, 
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Figure A. 3: These graphs show the NLTE departure coefficients bi on a logarithmic 
scale for all levels of oxygen against the layer number in the atmosphere model. The layer 
number increases going inwards. The number of levels in each ionization stage are written 
in parentheses. 

The wavelength variant (top) of the bound-bound rates and opacities, equations (A. 23) 
and (A. 22), are problematic in the inner regions of the atmosphere. There, the departure 
coefficients do not converge to bi = 1. That is the LTE approximation that holds better 
the deeper the layer in the atmosphere is. 

The frequency variant, equations (A. 15) and (A. 14), does not show this problematic be- 
havior and for this and other reasons, see the text, is superior to the other variant. 
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instead of assuming that all interactions of the radiation with the matter to- 
gether reproduce the TE radiation field, each transition on itself reproduces 
the TE field. This requirement is stronger than the requirement of detailed 
radiative balance. 

The exponential terms in (A. 25) originate from using the Boltzmann 
law, which is exactly the difference to equation (A. 10). Note that both 
and are strongly wavelength dependent. In the line center this form of 
the Einstein relations yields the wavelength variant of equations equation 
(A.6) and (A.7). 

As constant Aji coefficients are the line strength input data to PHOENIX 
these must be converted to B^. In order to do this the Aji are interpreted as 
the line center A°- values and converted to B}- using the Einstein relations 
for the line centers. 

In this variant only the Bf- are constant and therefore these are the 
appropriate coefficients to express the rates and opacities in 

R% = B^\ <fi\J\ dA (A.26) 
Jo 



RO. = B ^ e hc/XokT J°° + j^j ^-hc/XkT dA (A 27) 

^ = 4^^ n ^¥| e " CAfeT (A - 28) 



he 
47rA 



(A.29) 



4< (<M* " i>xb*e- hc ' XkT ) 



In this variant the wavelength- Voigt profiles are used, equations (A. 17) to 
(A. 19), as originally in PHOENIX. 

The single line source function for this original variant becomes 

bb;0 _ 2hc 2 ( b* <p X hc/XkT 1 \ rA o m 

When the Boltzmann distribution holds and complete redistribution is as- 
sumed then this reduces to the Planck function. 

The results with this variant for the same layer (deep in the atmosphere) 
as in figures A.l and A. 2 are shown in figure A. 4. The opacities (emission 
and extinction) are dominated by the bound-bound opacities shortwards of 
lOOA. This is partly due to the wavelength- Voigt profiles. Like discussed on 
page 149 these profiles exhibit a strong blue line wing. The mean intensity 
J\ accurately reproduces the Planck function. However, the accuracy is so 
high that in the second moment of the radiation field the accurate intensities 
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Figure A. 4: This figure shows the same quantities as figures A.l and A. 2 with the only 
difference that the bound-bound opacities are computed in the original PHOENIX variant, 
equations (A.29) and (A. 28). 

It produces an almost perfect Planckian source function for the whole wavelength range, 
because even the single line source functions, equation (A. 30), are Planckian as b* « 1. 
The integrated Eddington flux H is lower by a factor of 12 compared to the frequency 
variant. The reason is that the almost perfect Planckian intensities cancel efficiently in 
the second moment of the radiation field. 

In this variant wavelength- Voigt profiles are used, which feature strong blue line wings. 
These blue wings dominate the short wavelength range of the extinction coefficient. Due 
to the strong wavelength dependence of Aji the bound-bound emissivity looks Planckian. 



156 APPENDIX A. PROBLEMS WITH BROAD LINE THEORY 



efficiently cancel. For this original PHOENIX variant the integrated Eddington 
flux H for this layer is a factor of 12 smaller than for the frequency variant. 
On itself this is not yet a big problem. 

Big problems with this variant arise in the outer parts of the atmosphere. 
Figure A. 5 shows the same model as in figure A. 4 for a outer layer of the 
atmosphere. When for a radiative transition the departure coefficient of 
the upper level is larger than of the lower level b* > b* and the wavelength 
increases, then at some wavelength the extinction, equation (A. 29), becomes 
negative. Negative extinction on itself is not unphysical, this phenomenon 
is called lasering, see [Rut95]. However, lasering effects can cause the (to- 
tal) source function to become negative. A negative source function in the 
radiative transport formalism can lead to negative intensities along a char- 
acteristic ray. This is clearly unphysical. And numerically the operator 
splitting method does not converge. Therefore, if for a transition lasering 
would occur, the opacities of the transition are omitted. The extinction 
goes to zero smoothly with wavelength, but the emissivity does not at the 
same time. Omitting the transition longwards of the wavelength from which 
lasering occurs leads to jumps in the bound-bound emissivity. These jumps 
are the side effect of not treating lasering in the radiation transport. The 
physically right way to deal with this problem would be to find a way to 
treat the lasering effect in the radiative transport. This has not yet been 
accomplished in the scope of this work, but is suspended for future work. 

Alternatively, one could omit all opacity of lines with b* > b*. However, 
this condition is much stronger than the condition nj/gj > ni/gi that is used 
in the frequency variant. In LTE higher levels are energetically unfavored, 
described by the exponential damping term in the Boltzmann law, equation 
(2.12). Experiment shows that omitting such lines throws out a lot of lines 
that are treated in the Frequency variant, for example the vast majority of 
lines below 20 A. 

The frequency variant in the outer layers, figure A. 6, shows the same 
artifact at small wavelengths as in inner layers, figure A. 2. The emissivity 
does not fall off steeply in the Wien tail. In the frequency variant the only 
wavelength dependence in the emissivity comes from the profile function. 
However, the overall effect of this problem on the model is small. The 
rates are not affected significantly, because of the A 2 /c term being small for 
short wavelengths (equations (2.33) and (2.45)). And also the influence on 
the opacity integrals of the temperature correction (see section 2.11.1) is 
negligible. 
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Figure A. 5: This figure shows the same model as figure A. 4 for an outer layer in the 
atmosphere (layer 20) on a logarithmic wavelength scale. The bound-bound opacities are 
computed in the original PHOENIX variant, equations (A. 29) and (A. 28). 
Jumps in the emissivity occur due to omitting transitions that would exhibit lasering. 
Presently, lasering cannot yet be treated in the radiative transport. At the jump points 
the line extinction coefficient smoothly goes to zero, so the run of the extinction is smooth. 
The artifacts in the emissivity also leave their signatures in the flux. Comparison of the 
flux in this layer with the frequency variant, figure A. 6, shows that the frequency variant 
is preferable. 
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Figure A. 6: This figure shows the same quantities as figure A. 5 for the same atmo- 
spheric iayer on a iogarithmic waveiength scaie. The models are identical except that here 
the bound-bound opacities are computed in the frequency variant, equations (A. 14) and 
(A.15). 

The main differences lie in the shape of the emissivity and the flux. One disadvantage of 
the frequency variant is that due to unrealistic emission line profiles the emissivity does 
not fall of steep enough in the short wavelength range. This causes the mean intensity and 
the flux to rise again towards shorter wavelengths, which is clearly unphysical. However, 
the influence of this problem on the model is negligible (see the text). 



Appendix B 

First non-solar fits to V4743 
Sgr 



The models presented here are first-go, non-solar abundance fits to the five 
V4743 Sgr grating spectra that were also shown in section 5.4. Their com- 
putation just finished the very same day on which this work was finalized. 
Today, September 7th 2009, the new supercomputer in the Konrad-Zuse- 
Zentrum in Berlin became available and was yet completely empty, so that 
7x480 processors could be reserved for a few hours even before the opening 
was officially announced. With this total of 3360 processors, the last stage 
(stage 3, see section 4.3) of a whole model grid with Voo = 2400 km/s could 
be computed at once. 

The chosen abundances are solar, except for N and O, which are over- 
abundant by a factor of 100 and 10 respectively. These abundances are a 
first endeavor for non-solar abundances, based on the compositions obtained 
from nova evolution models for the ejecta [SIH + 09]. 

As yet, these abundances N lOOx solar and O lOx solar have not yet been 
tuned in any way. It is the first and only non-solar composition that was 
computed. Nevertheless, the models show a very interesting improvement 
of the fits to the observations with respect to the solar abundance models 
shown in section 5.4. 

Figures B.l to B.5 show the following observations. 
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Just like for the solar abundance models the terminal velocity is Voo = 2400 
km/s, and f3 = 1.5. The the non-solar models are plotted in green along 
with the solar models, for comparison, in red. 

The most significant improvements from the non-solar abundances with 
lOOx N and lOx O relative to solar, is the stronger N VII ionization edge at 
22.5AA and absorption lines at 24.8AA (Nvn) and 21.6AA (Ovu). 
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Figure B.l: Expanding model 




Figure B.2: Expanding model 
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Figure B.3: Expanding model 
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Figure B.4: Expanding model 




Figure B.5: Expanding model 



APPENDIX B. FIRST NON-SOLAR FITS TO V4743 SGR 



Bibliography 



[AGS+05] M. Asplund, N. Grevesse, A. J. Sauval, C. Allende Prieto, and 
D. Kiselman. Line formation in solar granulation. IV. [O I], 
O I and OH lines and the photospheric O abundance. A&A, 
435:339-340, May 2005. 

[Aue71] L. H. Auer. The stellar atmospheres problem. Journal of Quan- 
titative Spectroscopy and Radiative Transfer, 11:573-587, 1971. 

[AufOO] J. P. Aufdenberg. Line-blanketed spherically extended model 
atmospheres of hot luminous stars with and without winds. 
Ph.D. Thesis, August 2000. 

[Baa92] R. Baade. Line formation in multi-dimensional atmospheres with 
arbitary velocity fields. In C. S. Jeffery and R. E. M. Griffin, 
editors, Stellar Chromospheres, Coronae and Winds, pages 49- 
+, 1992. 

[BE08] M. F. Bode and A. Evans. Classical Novae. April 2008. 

[BKR+96] R. Baade, T. Kirsch, D. Reimers, F. Toussaint, P. D. Bennett, 
A. Brown, and G. M. Harper. The Wind Outflow of zeta Aurigae: 
A Model Revision Using Hubble Space Telescope Spectra. ApJ, 
466:979-+, August 1996. 

[BM92] M. Balucinska-Church and D. McCammon. Photoelectric ab- 
sorption cross sections with variable abundances. ApJ, 400:699- 
+, December 1992. 

[Boe89] E. Boehm-Vitense. Introduction to stellar astrophysics. Vol. 2. 
Stellar atmospheres. 1989. 

[CAK75] J. I. Castor, D. C. Abbott, and R. I. Klein. Radiation-driven 
winds in Of stars. ApJ, 195:157-174, January 1975. 

[Dir47] P. A. M. Dirac. The principles of quantum mechanics. 1947. 



167 



168 



BIBLIOGRAPHY 



[DLM+97] K. P. Dere, E. Landi, H. E. Mason, B. C. Monsignori Fossi, and 
P. R. Young. CHIANTI - an atomic database for emission lines. 
A&AS, 125:149-173, October 1997. 

[DLYD01] K. P. Dere, E. Landi, P. R. Young, and G. Del Zanna. CHIANTI- 
An Atomic Database for Emission Lines. IV. Extension to X-Ray 
Wavelengths. AJS, 134:331-354, June 2001. 

[Dre03] S. Dreizler. Temperature Correction Schemes. In I. Hubeny, 
D. Mihalas, and K. Werner, editors, Stellar Atmosphere Model- 
ing, volume 288 of Astronomical Society of the Pacific Confer- 
ence Series, pages 69 — h, January 2003. 

[Gri74] H. R. Griem. Spectral line broadening by plasmas. 1974. 

[Gro37] W. Grotrian. Zur physikalischen Deutung der Lichtkurve der 
Nova Herculis 1934. Mit 3 Abbildungen. Zeitschrift fur Astro- 
physik, 13:215-+, 1937. 

[GS78] J. S. Gallagher and S. Starrfield. Theory and observations of 
classical novae. ARA&A, 16:171-214, 1978. 

[Hau92] P. H. Hauschildt. A fast operator perturbation method for the 
solution of the special relativistic equation of radiative transfer 
in spherical symmetry. JQSRT, 47:433, 1992. 

[HB99] P. H. Hauschildt and E. Baron. Numerical solution of the ex- 
panding stellar atmosphere problem. Journal of Computational 
and Applied Mathematics, 109:41-63, September 1999. 

[HB04] P. H. Hauschildt and E. Baron. Improved discretization of the 
wavelength derivative term in CMF operator splitting numerical 
radiative transfer. A&A, 417:317-324, April 2004. 

[HB08] P. H. Hauschildt and E. Baron. A 3D radiative transfer frame- 
work. III. Periodic boundary conditions. A&A, 490:873-877, 
November 2008. 

[HBBA03] P. H. Hauschildt, T. S. Barman, E. Baron, and F. Allard. Tem- 
perature Correction Methods. In I. Hubeny, D. Mihalas, and 
K. Werner, editors, ASP Conf. Ser. 288: Stellar Atmosphere 
Modeling, pages 227 — h, January 2003. 



[HG03] W.-R. Hamann and G. Grafener. A temperature correc- 
tion method for expanding atmospheres. A&A, 410:993-1000, 
November 2003. 



BIBLIOGRAPHY 



169 



[HLB01] Yozo Hida, Xiaoye S. Li, and David H. Bailey. Algorithms for 
quad-double precision floating point arithmetic. In IEEE Sym- 
posium on Computer Arithmetic, pages 155-162, 2001. 

[HSB94] P. H. Hauschildt, H. Storzer, and E. Baron. Convergence prop- 
erties of the accelerated A-iteration method for the solution 
of radaitive transfer problems. Journal of Quantitative Spec- 
troscopy and Radiative Transfer, 51:875-891, June 1994. 

[HSS+95] P. H. Hauschildt, S. Starrfield, S. N. Shore, F. Allard, and 
E. Baron. The Physics of Early Nova Spectra. ApJ, 447:829-+, 
July 1995. 

[ISK + 00] N. Itoh, T. Sakamoto, S. Kusano, S. Nozawa, and Y. Kohyama. 

Relativistic Thermal Bremsstrahlung Gaunt Factor for the Intra- 
cluster Plasma. II. Analytic Fitting Formulae. AJS, 128:125-138, 
May 2000. 

[KPP99] J. Kubat, J. Puis, and A. W. A. Pauldrach. Thermal balance 
of electrons in calculations of model stellar atmospheres. A&A, 
341:587-594, January 1999. 

[KPPA89] R. P. Kudritzki, A. Pauldrach, J. Puis, and D. C. Abbott. 

Radiation-driven winds of hot stars. VI - Analytical solutions 
for wind models including the finite cone angle effect. A&A, 
219:205-218, July 1989. 

[Lam97] H. J. G. L. M. Lamers. The Theory of Line Driven Stellar Winds. 

LNP Vol. 497: Stellar Atmospheres: Theory and Observations, 
497:159-+, 1997. 

[LDY+06] E. Landi, G. Del Zanna, P. R. Young, K. P. Dere, H. E. Mason, 
and M. Landini. CHIANTI-An Atomic Database for Emission 
Lines. VII. New Data for X-Rays and Other Improvements. AJS, 
162:261-280, January 2006. 

[Luc64] L. B. Lucy. A Temperature-Correction Procedure. SAO Special 
Report, 167:93-+, December 1964. 

[MC37] D. H. Menzel and G. G. Cillie. Hydrogen Emission in the Chro- 
mosphere. ApJ, 85:88-+, March 1937. 

[McL57] D.B. McLaughlin. Stars and Stellar Systems. 1957. 

[Mih70] D. Mihalas. Stellar atmospheres. Series of Books in Astronomy 
and Astrophysics, San Francisco: Freeman, — cl970, 1970. 



[Mih78] D. Mihalas. Stellar atmospheres /2nd edition/. 1978. 



170 



BIBLIOGRAPHY 



[Mih80] D. Mihalas. Solution of the comoving-frame equation of trans- 
fer in spherically symmetric flows. VI - Relativistic flows. ApJ, 
237:574-589, April 1980. 

[MW84] D. Mihalas and B. Weibel Mihalas. Foundations of radiation 
hydrodynamics. 1984. 

[Nes09] J.-U. Ness, privat communication, 2009. 

[NSB+03] J.-U. Ness, S. Starrfield, V. Burwitz, R. Wichmann, 
P. Hauschildt, J. J. Drake, R. M. Wagner, H. E. Bond, J. Kraut- 
ter, M. Orio, M. Hernanz, R. D. Gehrz, C. E. Woodward, 
Y. Butt, K. Mukai, S. Balman, and J. W. Truran. A Chandra 
Low Energy Transmission Grating Spectrometer Observation of 
V4743 Sagittarii: A Supersoft X-Ray Source and a Violently 
Variable Light Curve. AJL, 594:L127-L130, September 2003. 

[NSB+07] J.-U. Ness, S. Starrfield, A. P. Beardmore, M. F. Bode, J. J. 

Drake, A. Evans, R. D. Gehrz, M. R. Goad, R. Gonzalez-Riestra, 
P. Hauschildt, J. Krautter, T. J. O'Brien, J. P. Osborne, K. L. 
Page, R. A. Schdnrich, and C. E. Woodward. The SSS Phase 
of RS Ophiuchi Observed with Chandra and XMM-Newton. I. 
Data and Preliminary Modeling. ApJ, 665:1334-1348, August 
2007. 

[OK87] G. L. Olson and P. B. Kunasz. Short characteristic solution of 
the non-lte line transfer problem by operator perturbation, i. - 
the one-dimensional planar slab. JQSRT, 38:325, 1987. 

[OPG+06] J. Osborne, K. Page, A. B. M. Goad, M. Bode, T. O'Brien, 
G. Schwarz, S. Starrfield, J.-U. Ness, J. Krautter, J. Drake, 
A. Evans, and S. P. S. Eyres. RS Oph: Swift X-ray observations 
find short period modulation and highly variable low energy flux. 
The Astronomer's Telegram, 770:1 — h, March 2006. 

[Pay57] OH. Payne. The Galactic Novae. 1957. 

[Pet05] A. Petz. Modeling atmospheres of classical novae in X-rays with 
PHOENIX. PhD thesis, Hamburger Sternwarte, Oktober 2005. 

[PHNS05] A. Petz, P. H. Hauschildt, J.-U. Ness, and S. Starrfield. Mod- 
eling CHANDRA low energy transmission grating spectrometer 
observations of classical novae with PHOENIX. I. V4743 Sagit- 
tarii. A&A, 431:321-328, February 2005. 



[Pom73] G. C. Pomraning. The equations of radiation hydrodynamics. 
1973. 



BIBLIOGRAPHY 



171 



[PSS78] D. Prialnik, M. M. Shara, and G. Shaviv. The evolution of a 
slow nova model with a Z = 0.03 envelope from pre-explosion to 
extinction. A&A, 62:339-348, January 1978. 

[PSS79] D. Prialnik, M. M. Shara, and G. Shaviv. The evolution of a 
fast nova model with a Z = 0.03 envelope from pre-explosion to 
decline. A&A, 72:192-203, February 1979. 

[RL79] G. B. Rybicki and A. P. Lightman. Radiative processes in astro- 
physics. New York, Wiley-Inter science, 1979. 393 p., 1979. 

[Rut95] R. Rutten. Radiative transfer in stellar atmospheres. 1995. 

[SBLR01] R. K. Smith, N. S. Brickhouse, D. A. Liedahl, and J. C. Ray- 
mond. Standard Formats for Atomic Data: the APED. In 
G. Ferland and D. W. Savin, editors, Spectroscopic Challenges 
of Photoionized Plasmas, volume 247 of Astronomical Society of 
the Pacific Conference Series, pages 161 — h, 2001. 

[SIH+09] S. Starrfield, C. Iliadis, W. R. Hix, F. X. Timmes, and W. M. 

Sparks. The Effects of the pep Nuclear Reaction and Other Im- 
provements in the Nuclear Reaction Rate Library on Simulations 
of the Classical Nova Outburst. ApJ, 692:1532-1542, February 
2009. 

[SS87] S. Starrfield and W. M. Sparks. A Review of the Classical Nova 
Outburst. APSS, 131:379-393, March 1987. 

[SSMS97] M. Steffen, R. Szczerba, A. Men'shchikov, and D. Schoenberner. 

Hydrodynamical models and synthetic spectra of circumstellar 
dust shells around AGB stars. A&AS, 126:39-65, November 
1997. 

[SST74] S. Starrfield, W. M. Sparks, and J. W. Truran. CNO abundances 
and hydrodynamic models of the nova outburst. II - 1.00 solar 
mass models with enhanced carbon and oxygen. AJS, 28:247- 
270, August 1974. 

[Sta89] S. Starrfield. Thermonuclear processes and the classical nova 
outburst. In Classical Novae, pages 39-60, 1989. 

[STWS96] S. Starrfield, J. W. Truran, M. Wiescher, and W. M. Sparks. Nu- 
cleosynthesis and the Nova Outburst. In S. S. Holt and G. Son- 
neborn, editors, Cosmic Abundances, volume 99 of Astronomical 
Society of the Pacific Conference Series, pages 242 — h, 1996. 

[Sut98] R. S. Sutherland. Accurate free- free Gaunt factors for astrophys- 
ical plasmas. MNRAS, 300:321-330, October 1998. 



172 

[Und49] 
[Uns38] 
[Uns55] 
[WAMOO] 

[WZ72] 



BIBLIOGRAPHY 



A. B. Underhill. On the Effect of Radiation Pressure in the 
Atmospheres of Early-Type Stars. MNRAS, 109:562-+, 1949. 

A. Unsold. Physik der sternamosphdren, MIT besonderer beruck- 
sichtigung der Sonne. 1938. 

A. Unsold. Physik der Sternatmosphdren, MIT besonderer 
Beriicksichtigung der Sonne. 1955. 

J. Wilms, A. Allen, and R. McCray. On the Absorption of X- 
Rays in the Interstellar Medium. ApJ, 542:914-924, October 
2000. 

J. W. Wijbenga and C. Zwaan. Empirical NLTE Analyses of 
Solar Spectral Lines. I: A Method and Some Applications to 
Earlier Analyses. Sol.Phys., 23:265-286, April 1972. 



Acknowledgements 



First, I am delighted to thank Prof. Peter Hauschildt for providing me with 
the very interesting topic, for the lots of good advice, especially in the last 
months, and for helping me with the successfull execution of this work. 

I thank Prof. Eddie Baron for the very successful collaborations that we 
had. I look forward to meet you in Chicago. 

I thank Prof. Sumner Starrfield for the short and inspiring talk we had 
in Tempe. There will be a lot of good nova work to do together. 

A warm thanks to Knop, my very nice and very helpful roommate, I 
thank you for all the help you gave me, for taking the time to discuss prob- 
lems with me. 

Thanks to Jan-Uwe Ness, for lots of good collaboration, for providing me 
with observational data, for your hospitality, for the interesting insights into 
the scientific world and the nova community, for introducing me to people, 
motivating me to attend workshops, write papers etc. 

Also I am very thankful to Prof. Rob Rutten for the incredibly inspiring 
graduate course and exercises in radiative transport that he gave in Utrecht 
2001. While this was the most interesting course I ever had on university, 
he gave me the fundamental understanding that I'm still profited by. 

I thank Dr. John Dombeck for his great idl color tables. 

Then, I want to thank Oom Piet and Tante Cora for the many encour- 
aging talks, for all their support for us, the good and useful advice out of 
long good experience in a situation like ours. 

I thank my parents for all their very important advice for important 
decisions that were to be made lately, and in the past. 

Finally, the person that I want to thank most is my lovely wife Silke, for 
supporting me with endless patience in the lots of work that was to be done 
in the last couple of years, months, weeks and days! 



173 



Erklaerung 



Hiermit versichere ich, Daniel R. van Rossum, die vorliegende Arbeit selb- 
staendig verfasst und nur unter Zuhilfenahme der angegebenen Quellen und 
Hilfsmittel angefertigt zu haben. Desweiteren erklaere ich mich mit dem 
Verleih und der Veroeffentlichung dieser Arbeit einverstanden. 



Unterschrift 



Hamburg, den 31.05.2006 



